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ruination  of  mean  gravity  anomalies  at  the  surface  of  the  earth  using  least 
squares  collocation  and  theoretical  covariance  functions.  This  paper 
discusses  the  theoretical  basis  of  the  above  procedures  and  describes  two 
simulation  studies  performed. 

*KS 

In  the  first,  postulated  data  surrounding  a 5°  equal  area  block  and  a 
2°  x 2°  block  was  used  to  estimate  the  accuracy  in  which  the  anomaly  could 
be  recovered.  For  example,  with  acceleration  data  having  an  accuracy  of 
±1  mgal,  a 5°  anomaly  could  be  determined  with  an  accuracy  of  ±4  mgals 
with  the  low  satellite  at  250  km,  and  ±8  mgals  with  the  low  satellite  at 
850  km.  A comprehensive  simulation  experiment  was  then  performed  to 
check  the  actual  recovery  of  postulated  anomalies  determined  from  defined 
sets  of  potential  coefficients.  Assuming  known  orbits  the  recovery  of  the 
unknown  anomalies  was  accomplished  to  about  ±2  mgals.  When  orbits 
errors  were  introduced  the  errors  increased  but  could  properly  be  con- 
trolled through  assignment  of  data  accuracies. 

The  promising  results  of  this  study  indicate  that  least  squares 
collocation  techniques  can  be  advantageously  applied  to  this  type  of  anomaly 
recovery  avoiding  the  instability  of  the  downward  continuation  problem  that 
exists  in  other  methods  of  anomaly  recovery  from  this  data  type. 
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0.  Introduction 

Our  goal  is  the  derivation  of  gravity  anomalies  at  the  surface  of  the  earth 
from  satellite  to  satellite  range  rates.  The  determination  of  gravity  parameters 
from  SST-observations  is  described  in  several  reports,  among  others  in 
(Schwarz,  1970;  Kaula,  1972;  Martin,  1972;  Hajela,  1974a;  Kaula  et  al,  1975  ). 
Common  to  the  techniques  described  there  is  that  the  partial  derivative  of  the 
sum  range  and  range  rate  are  formed  with  respect  to  the  unknown  parameters. 
From  the  partials  a design  matrix  is  built  and  the  residual  ranges  and  range 
rates  (observed  minus  computed)  are  modeled  in  the  least  squares  sense  by 
the  unknown  parameters.  Since  implicitly  in  all  methods  of  surface  gravity 
parameter  determination  from  satellite  observations  a downward  continuation 
problem  is  involved,  certain  difficulties  arise  in  the  mentioned  least  squares 
methods.  These  difficulties  are  in  keywords:  truncation  of  the  parameters 
to  a "inner  cap",  stability  of  the  solution,  definition  of  resolution  and  so  on. 
Therefore,  Rummel  (1975a,  1975b)  proposed  the  method  of  least  squares 
collocation  for  the  stated  problem.  Least  squares  collocation  should  avoid 
the  mentioned  problems  to  a great  extent.  Since  it  will  be  split  in  different 
steps,  each  step  can  be  separately  checked  for  certain  criteria  such  as 
stability  and  resolution. 

In  the  sequel  the  method  of  determination  of  gravity  anomalies  by  least 
squares  collocation  is  described.  The  method  will  in  a certain  sense  combine 
some  ideas  of  Muller  and  Sjogren  (1968),  and  of  Kaula  (1972)  with  least 
squares  collocation. 

In  this  report  a theoretical  procedure  for  recovery  of  gravity  anomalies 
from  SST  data  will  be  given  in  detail.  The  theory  will  then  be  examined 
through  simulation  studies  that  will  demonstrate  the  accuracy  of  the  pro- 
posed method  and  indicate  the  accuracy  to  be  expected  of  the  recovered 
anomalies. 


1.  Theoretical  Discussion 

We  are  considering  satellite  to  satellite  tracking  (SST)  observations 
in  a "high-low"  experiment  such  as  between  the  ATS-6  and  GEOS-3 
satellites.  The  data  gathered  in  a SST  experiment  are  sum  range  and 
sum  range  rate,  both  types  being  described  in  (Martin,  1972;  Hajela,  1974). 
Because  of  the  low  information  content  of  sum  range  data  for  gravity 
parameter  determination  only  sum  range  rate  data  will  be  analyzed  in 
this  study. 
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The  observed  sum  range  rate,  Rob< , 


is  defined  by: 


(1.1) 


Rot  • 


+ Rrtd 


+ Rr  c u Rfi  r i 


with  subscripts  g,  c,  and  r for  ground  tracking  station,  close  satellite 
(e.g.  G EOS-3),  and  relay  satellite  (e.g.  ATS-6),  and  with  subscripts 
u and  d for  measurements  directed  upwards  and  downwards. 

Because  the  gravity  field  in  the  altitude  of  the  relay  satellite  is 
highly  attenuated,  the  range  rates  ru  and  R5r4  may  be  computed 
very  well  from  a low  degree  set  of  potential  coefficients.  Doing  this, 
we  find  the  mean  range  rate  ftr0  between  the  relay  and  the  close 
satellite  from: 

(1.2)  Rrc  = \ ( Rr<*i  + Rrou)  = Rob.  -i(R*ru  +R*rd/) 


The  time  variation  of  ftre  is  expressed  by: 


(1.3) 


d Rr  o 

dt 


ARro 

At 


The  acceleration  Rr  0 is  central  for  our  computations  and  will  now  be 
discussed  in  more  detail. 


The  orbit  of  the  relay  as  well  as  of  the  close  satellite  is  defined  at 
every  instant  by  their  position  vectors  X and  velocity  vectors  X-  After 
applying  corrections  for  side  effects  such  as  solar  radiation  pressure 
and  atmospheric  drag  the  acceleration  vector  X is: 


(1.4)  X = VV  , 

where  W is  the  gradient  of  the  actual  harmonic  potential  of  the  earth. 
The  difference  of  the  acceleration  vector  of  the  relay  satellite,  and 
that  of  the  close  satellite,  X_= » yields: 

(1.5)  Xro  - _ ^Vo  % 
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The  doppler  instrumentation  detects  only  the  quantity  f£  0 of  equation 
(1.3),  which  is  the  magnitude  of  the  vector  Rr  o in  a straight  line  connecting 
the  locations  S,  and  S,  of  the  two  satellites.  In  other  words,  RrC  represents 
the  projection  of  the  vector  Xj  * onto  the  line  ST^e  • Denoting  the  unit 
vector  directing  from  Sr  to  S0  by  » we  obtain: 

Rr0  = Xrc  * eye,  or  using  equation  (1.5) 

(1.6)  kr o = vVr  * e^o  ~ ¥¥.<>  ' 

These  basics  are  demonstrated  in  Figure  1.1. 

As  already  mentioned,  the  motion  of  the  relay  satellite  can  be  well 
described  from  a low  degree  set  of  potential  coefficients.  That  means  that 
the  gradient  VVr  can  be  determined  from  this  set.  Computing  in  addition 
a reference  gradient  VVT'for  the  close  satellite  a reference  acceleration 
h*?of  is  obtained  from: 

rT?  = tvr  - ^yrf)  * 

With  the  anomalous  potential,  T,  defined  by  T=V-Vr9f  we  derive  the  residual 
acceleration,  RjV  • from: 


Equations  (1. 1)  to  (1. 7)  provide  the  relation  of  the  observed  range 
rates  with  the  gradient  of  anomalous  potential.  The  gravity  anomalies, 

Ag,  on  the  surface  of  the  earth  may  now  be  estimated  from  the  gradient 
of  the  anomalous  potential  at  the  altitude  of  the  close  satellite  by  a 
straightforward  application  of  least  squares  collocation.  Only  the  fact 
that  the  input  quantity  to  collocation  is  in  this  case  a vector  process  might 
deserve  some  additional  thoughts.  The  least  squa-os  collocation  estimate 
of  Ag  at  any  point  Q on  the  surface  of  the  earth  is  (Rummel,  1975a): 
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p.l 


(1.8)  tfg  (Q)  = C< Q)AtV!|  (C  + Djif1  Zlt 

dll  1*1  l»3a  3»*3n  3 n • 3 n 3n*l 

with  C(Q)^8^ij  ....  crosscovariance  vector  between 

Ag(Q)  and  the  input  elements  7T , 


at  the  points  P,  j=i,  ...n,  dim(l-3n), 

C 2l  j 21  matrix  of  the  autocovariances  between 

the  anomalous  potential  gradients 
2T3  with  dim(3q*3n), 

D } i noise  matrix 


What  we  would  need  to  apply  this  equation  to  anomaly  recovery  is  the  vector 
of  gradients  vTt . But  actually  we  have  only  fVc®6 » the  projection  of  the 
gradient  onto  the  line  connecting  the  close  and  the  relay  satellite  as  expressed 
in  equation  (1.7).  In  other  words  we  would  like  to  find  the  inverse  of  this 
projection,  that  means  an  expression  of  the  form: 

VT0  = Afi”e‘ 

In  order  to  define  the  projector  A one  would  have  to  know  not  only  the 
direction  defined  by  ere  but  also  the  direction  of  the  gradient  vTt  , with 
unit  vector,  e yj.  , which  is  not  available  to  us.  The  problem  becomes 
obvious  by  going  back  to  equation  (1.  7):  Rewritting  equation  (1.  7)  in 
components  of  a rotating  (X2  , Xg,  X3  ) coordinate  system  we  obtain 


5Tr.  Xi 
3Xl  - 


Sfc-*"' 


^ rr 

This  equation  contains  the  three  unknowns,  ° * i 33,  which  means  that 
we  need  at  least  three  linearly  independent  1 equations  of  this  type  at 
every  instant  in  order  to  be  able  to  determine  the  gradient  VTC . As  a result 
we  note: 


For  a unique  recovery  of  gravity  quantities  from  a"ligh-low"  SST 
experiment  at  least  three  relay  satellites  are  needed  that  track  the 
close  satellite  simultaneously. 


Instead  of  postulating  this  ideal  situation  we  will  assume  the  more  realistic 
case  of  only  a satellite  pair.  From  the  above  result  it  becomes  obvious 
that  a certain  approximation  has  to  be  introduced  at  this  point.  Basically  we 
may  choose  between  two  types  of  approximations,  both  with  respect  to  the 
direction  of  the  gradient,  VT0  . 
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Approximation  (I): 

(1.9)  VTe  = ( , 0 , o)T  ((.  )T  ...  transposed  ) 

' ^ r 

which  means  the  gradient  of  the  anomalous  potential  has  approximately  only  a 
radial  component. 

Approximation  (II): 

(1.10)  VTe  = vVT  , 


which  means  the  gradient  ?Te  has  approximately  the  direction  of  The 

geometrical  interpretation  for  both  approximations  is  given  in  Fig.  1.2. 

We  expect  that  by  using  approximation  (13)  a smaller  error  is  committed 
than  by  using  approximation  (I).  But  this  advantage  has  to  be  paid  by  a far 
higher  computational  effort  since  (1)  for  each  observational  point  TVS'*  has 
to  be  computed  and  (2)  multivariate  collocation  has  to  be  used  as  expressed 
by  equation  (1.8). 

Since  we  are  convinced  that  approximation  (I)  is  for  most  practical  situ- 
ations suf  ient  we  continue  to  analyze  only  the  latter.  Assuming  equation 
(1.9)  to  be  valid  the  relation  to  the  acceleration  Tt! ? is: 


7T.  = .o.o  y = & 

— \ > r ' cc 


cos  8 


ec 


where  8 is  the  spatial  angle  between  the  vectors  VTC  and  fgV  with 

cos 8 = ec  • ero,  with  e0  the  unit  vector  in  the  radial  direction  at  the  close 

satellite. 

Using  only  the  radial  component  this  equation  becomes: 


(1.11) 


Tr 


U.  = J£L 

> r cos ft 


From 
to  errors 


equation  (1.11)  it  follows  that  errors,  eV  , in  r£®*  propagate 


•*t 

} r 


in 


vr_ 

a r 


by 


cos*  8 


var 


(1.12) 


™r  (to  ) 

*r 


j 


I 

! 


This  relation  demonstrates  the  importance  of  the  geometry  of  the  relay 
to  the  close  satellite  on  the  errors.  From  (1. 12)  follows; 

(1)  We  expect  good  results  from  an  analysis  of  acceleration  data 

at  times  where  the  close  satellite  was  moving  in  the  subsatellite 
region  of  the  relay  satellite. 

(2)  We  would  like  to  see  the  estimation  of  gravity  quantities  from 
SST-data  to  be  an  almost  local  procedure.  This  means  it  would 
be  desirable  that  only  data  gathered  in  the  subsatellite  region 
of  the  relay  satellite  would  have  to  be  introduced  into  the 
estimation  p roc  es  s . 

Equation  (1.11)  is  the  desired  relation  between  Rj“  as  obtained  from 
SST  and  dT  which  can  be  introduced  into  the  collocation  solution. 

Inserting  (1.11)  into  equation  (1.8)  yields: 


(1.13)  Ag(Q)  = £ (Q)^,,  (£Tr3Trl  + DjO'NRrT/cOSS): 


at*  i*i 


1 • n 


n*  n 


n *1 


Correspondingly  for  approximation  (II)  the  gradient  vTc  of  the  anomalous 
potential  is  related  to  Rj?  with  the  approximate  direction  of  vy0I‘8f  by; 


(1.14) 


VTe  = 


cos/3  * 


where  B is  now  the  spatial  angle  between  0 and  the  direction  of  vv£rf 
and  e v y is  the  unit  vector  of  . We  obtain  for  the  multivariate 
least  squares  collocation  equation  that  solves  for  Ag  by  inserting  (1.14) 
into  equation  (1. 8) 


(1.15) 


Ag(Q)  = C (Q)^g  j (£yvjyvi  + £li)x 


at*  1*1  1*  n 


3n  *3n 


a n*3n 


3t>'  1 


Computational  Procedure 

The  information  provided  to  the  user  in  a SST  experiment  is  usually 
the  time  of  measurement,  the  sum  range  and  sum  range  rate  observation, 
Robi  and  Rob,  respectively,  as  well  as  the  orbital  elements  at  certain 
periods  of  time.  Using  the  GEODYN  (Martin,  1972)  program  we  can 
compute  the  range  rate,  R^  9BP  , at  the  times  of  measurement  from  a set 
of  potential  coefficients  and  the  given  orbital  elements.  The  difference 
between  observed  and  computed  range  rate  is : 

(1.16)  f£V  = Rob,  - . 

which  is  the  relative  velocity  between  the  close  and  the  relay  satellite 
in  the  line  connecting  the  two  satellites.  As  indicated  in  equation  (1.3) 
we  have  then  to  find  the  acceleration  f£*e* . The  necessary  numerical 
differentiation  may  be  critical  for  the  whole  processing  and  therefore 
deserves  special  attention. 


We  propose  three  different  types  of  numerical  differentiation. 
Even  though  their  features  might  be  analyzed  theoretically,  only  their 
application  on  real  data  will  give  a final  answer  what  procedure  should 
be  chosen.  The  first  two  methods  of  numerical  differentiation  are 


deterministic.  This  is  to  a certain  extent  a drawback  because  the 
range  rate  data  will  be  disturbed  by  observational  errors.  But  after 
filtering  the  raw  data  these  procedures  should alsobe  applicable  with 
advantage. 


(1)  Finite  Differences;  (e.g.  Whittaker  and  Robinson,  1967,  p.  62) 
It  is  assumed  that  the  residual  range  rates,  r£*s  , are  given 
free  of  noise  equidistant  with  time  in  intervals  of  w.  Then, 
we  find 


T*' 


where  A,  As , A*,  A4  are  the  first  to  fourth  differences 
obtained  from  the  data  in  a difference  scheme.  Thus,  from 
each  sequence  of,  forexample  six,-- non  overlapping — data  an 
estimate  for  Rj*  can  be  determined. 


t 


(2)  Spline  Method:  (e.  g.  Shampine  and  Allen,  1973,  p.  54) 

This  method  seemed  more  accurate  to  us  than  method  (1)  and 
was  therefore  used  in  the  simulation. 

A cubic  interpolatory  spline  is  computed  for  each  sequence  of 
five  or  more — non  overlapping — R£7  data  by  a program  given 
in  (ibid).  By  a minor  modification  of  this  FORTRAN  program 
the  derivative  Rj?  is  obtained  at  any  desired  point  inside  the 
used  data  interval.  For  details,  see  Section  4.3, 

(3)  Least  Squares  Collocation: 

The  advantage  of  least  squares  collocation  would  be  that  in 
contrast  to  methods  (1)  and  (2)  the  filtering  of  observational 
noise  cau  be  implemented  into  the  procedure. 

Using  the  classical  least  squares  prediction  formula 

(1.17)  f£7  (t)  = C(t)tl  (Cttt3  + £t  tt  4 f1  f£o,  1 3 

we  obtain  from  the  vector  of  residual  range  rates  an  approxima- 
tion to  the  range  rate  function  Rj®0‘  (t)  at  any  desired  point  on  the 
real  (time-)  line.  The  quantities  appearing  in  equation  (1. 17)  are: 
C (t)t  j . . . the  autocovariance  vector  between  the  observations 
and  the  prediction  point, 

Ctlt4  . . .the  autocovariance  matrix  of  the  observations,  and 
Dt  lt  j . . .the  observational  noise  matrix. 

The  autocovariance  function  C(t)  of  Rjf  can  be  determined  by 
approximating  to  an  empirical  covariance  sequence — as  derived 
from  simulated  noise  free  data —a  positive  semidefinite  function 
such  as  a Gaussian,  second-order  Markov,  or  third-order 
Markov  model  shown  in  (Jordan,  1972).  The  autocovariance 
function  C(t)  + D(t)  can  be  derived  in  the  same  way  using  real 
data.  From  thesetwo  functions, C (t)  and  C(t)  + D(t)  the  co- 
variance  vector  and  matrix  needed  in  equation  (1. 17)  can  be 
built. 

The  acceleration,  , is  now  computed  at  any  desired  point  on 
the  real  line  by  differentiating  equation  (1.17), 

(1.18)  + «, 


Since  the  differentiation  is  directly  carried  out  on  the  analytical 
expression  for  C (t)  equation  (1. 18)  can  be  interpreted  as  an  analytical 
differentiation  of  the  approximate  function,  ^“'(t).  derived  from  (1.16). 

The  solution  proposed  and  applied  by  Muller  and  Sjogren  (1968)  for  the 
lunar  experiment  is  of  type  (2). 

After  computing  R£",  cos 5 has  also  to  be  determined.  For  approx- 
imation (I)  we  have  cosB=  ec  • e,.  e where  ej.  „ is  computed  from  the 
geocentric  coordinates  of  the  relay  and  the  close  satellite  and  ec  from 
the  geocentric  coordinates  of  only  the  close  satellite.  The  geocentric 
coordinates  are  provided  by  the  GEODYN  program  together  with  the 
computed  range  rates. 

More  tedious  is  the  computation  of  cosff  using  approximation  (II). 

Then  we  use  cos£=eyy  • e ye  and  the  unit  vector  e ^ for  the  gradient 
of  the  reference  potential,  vyer®f  has  to  be  determined  for  each  instant 
from  the  set  of  potential  coefficients. 

Finally,  least  squares  collocation  in  form  of  equations  (1.13)  or 
(1.15)  is  applied  on  the  approximate  values  for  VTa  derived  from 
equations  (1. 11)  or  (1. 14).  The  covariance  elements  will  be  computed 
from  global  covariance  expressions  such  as  those  derived  in  (Tscherning 
and  Rapp,  1974).  Thereby  we  use  the  FORTRAN  programs  developed  by 
Tbchering  (1976)  which  provide  the  covariance  of  the  anomalous  potential, 
the  gravity  anomaly  and  their  first  and  second  derivatives.  The  es- 
timated gravity  anomalies  are  chosen  to  be  located  on  a sphere  whose 
center  is  at  the  earth's  center  of  mass  and  with  the  mean  radius  of  the 
earth.  Mean  anomalies  are  computed  by  numerical  integration  of  the 
appropriate  covariances. 

The  problem  of  determining  gravity  anomalies  at  the  surface  of  the 
earth  from  * T values  in  satellite  altitude  lies  in  the  fact  that  a downward 
^r 

continuation  has  to  be  carried  out.  This  will  give  us  two  types  of  difficulties: 
(1)  Small  errors  at  satellite  altitude  will  blow  up  to  large  errors  in  Ag 
after  the  downward  continuation  operation  is  applied  on  them.  (2)  The 
downward  continuation  operation  may  cause  severe  stability  problems. 

To  get  an  idea  about  this  operation  the  expression  C(Q)^t  tp  . 

(C_fr  i Tr  j + D)-1  eQuati°n  (1. 13)  has  been  computed  for  a point,  and  a 
l!5®  equal  area  block  that  is  to  be  estimated.  The  above  covariance 
expression  transforms  the  (Rj*e*  /cosfl)-values  to  Ag-  values  at  the 
surface  of  the  earth.  The  given  data  are  assumed  to  have  a spherical 
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distance  of  0°  to  90°  in  5°  steps  from  the  center  of  the  surface  ^g- block 
and  are  located  in  250  km  altitude.  The  noise  is  considered  to  be  zero. 
The  results  are  given  in  Figure  1. 3 , where  the  values  of  a vector 
A = C (0)A,Tr  j Cfr'jr,.  t are  presented  with  increasing  spherical  distance 
from  the  estimation  point.  Even  though  we  know  that  theoretically  we 

would  have  to  have  given  the  — — values  as  a continuous  function  on  a 
sphere  concentric  with  the  cen^r  of  the  earth  and  in  satellite  altitude, 
Figure  1.3  shows  that  for  our  practical  situation  data  have  only  to  be 
given  up  to  a maximum  spherical  distance  of  about  20°  from  the 
estimation  point.  For  a spherical  distance  of  10°  the  influence  of  the 
data  is  already  below  20%  of  the  influence  of  the  data  with  0°  spherical 
distance.  This  question  will  be  analyzed  further  in  detail  in  the 
simulation  study.  In  the  simulation  only  approximation  (I)  will  be 
considered. 


2.  The  Covariance  Functions 

The  evaluation  of  (1. 13) requires,  in  part,  the  determination  of  the 
covariances  between  the  3T  values  at  points  i and  j,  and  the  covariance 

* r 5X 
between  the  anomaly  being  predicted  and  the  observed  or  given  — — 

o r 

value.  Such  covariance  functions  can  be  derived  in  an  analytical  way 
using  an  assumed  model  for  the  anomaly  degree  variances  (Tscherning 
and  Rapp,  1974,  Tscherning,  1976).  In  the  Tscherning (1976)  dis- 
cussion covariances  were  computed  for  gravity  anomalies,  -3T  / and 

/r 

other  quantities  not  of  specific  interest  for  this  report.  We  let; 


A,  =L  or  -rA 

r 3 r 3r 


Then  the  covariance  between  the values  is: 

3 r 

(2.2)  cov(^-  | |^  | ) = M{(-rP  A,)  (-r„  A,  )} 

*r  p o r q 


= rP  ijj  cov(AP  Aq) 

Here  M indicates  an  average  value  over  the  earth  and  r is  the  distance 
from  the  center  of  the  earth  to  the  point  in  question. 
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with  Tr  assumed  to  be  located  in  250  km  altitude 


I 


The  value  of  cov(AP  Aq)  is  available  from  subroutine  COVAX(Tschemingf 
1976)  in  Eotovos  units  squared.  IE3  is  equal  to  1 x lO^mgal2/ m2  so  that 
the  units  of  the  left  hand  side  of  equation  (2. 2)  will  be  in  mgal2  if  cov(ApAq)  is 
multiplied  by  10-8  and  the  r values  are  given  in  meters. 

Using  a similar  procedure  we  can  write  for  the  covariance  between  the 
anomaly  at  point  P and  the  radial  derivative  of  the  disturbing  potential  at  Q-. 

T 

(2.3)  cov(Agp  , — | ) = -rq  cov(  AgP  Aq ) 

^ r Q 

where  cov(  AgP  Aq)  is  given  in  units  of  mgal*  E units  from  COVAX.  If  r is 
given  in  meters  the  units  on  the  left  hand  side  of  equation (2. 3)  will  be  in  mgal2  if 
cov(Agp  Aq)  is  multiplied  by  10“4.  Typical  numerical  values  for  these 
covariance  functions  are  given  in  Table  2.1  (for  an  ellipsoid  reference 
field)  and  Table  2. 2 (for  a degree  12  reference  field). 

Table  2.1 

Sample  Covariances  Needed  in  SST  Reduction 
Given  With  Respect  to  an  Ellipsoid  Reference  Field 
(Units:  mgal3) 


/*T  1 
cov  l 

W |p 

ZL  1 ) 
* r |Q ' 

cov(Agp 

3T 

9r 

\ 

J 

hP  =250  km 

hp  =850  km 

hp  =0  m 

hp  = 

0 m 

0° 

h,  =250  km 

hj  =850  km 

hq  =250  km 

h?  = 

850  km 

0 

327 

98 

-341 

-128 

1 

324 

97 

-329 

-128 

2 

316 

97 

-300 

-125 

5 

276 

93 

-216 

-112 

7 

247 

89 

-176 

-101 

10 

208 

82 

-134 

-84 

13 

175 

73 

-103 

-68 

15 

154 

67 

-86 

-59 

20 

110 

51 

-54 

-39 
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Table  2.2 

Sample  Covariances  Needed  in  SST  Reduction 
Given  With  Respect  to  a Degree  12  Reference  Field 
(Units:  mgal8) 


Examination  of  the  values  given  in  Tables  2.1  and  2.2  show  the  rapid 
decrease  in  correlation  as  the  altitude  increases  from  250  km  to  850  km. 
This  decrease  is  simply  a numerical  representation  of  the  attentuation  of 
the  anomalous  gravity  field  with  height.  We  also  note  that  the  covariances 
with  respect  to  the  ellipsoidal  field  change  very  slowly  decreasing  to  the 
first  zero  crossing  for  the  radial  derivative  autocovariances  at  = 38° 

(for  both  altitudes)  and  for  the  covariance  between  the  surface  anomaly  and 
the  radial  derivative  at  i/i  = 34°  (for  both  altitudes).  When  the  higher  degree 
surface  is  used  as  reference  the  magnitude  of  the  covariances  considerably 
decreases  and  the  first  zero  crossings  occur  at  smaller  y values.  This 
would  indicate  that  only  data  close  to  the  computation  area  can  be  used 
when  the  higher  degree  reference  surface  is  used.  This  statement  will 
be  examined  in  the  simulation  studies  to  be  discussed  in  the  next  sections 
of  this  report. 


3.  Preliminary  Simulation  Studies 

Before  extensive  simulation  studies  are  performed  it  is  of  interest 
to  investigate  the  accuracy  to  be  expected  in  the  recovery  of  gravity 
anomalies  from  5T/9r  values.  Using  the  notation  of  equation  (1. 13) 
the  error  variance  of  the  estimated  gravity  anomaly  is  : 
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(3.1)  m2  ^Ag  (Q)  ) = C ( Ag ) - C(Q)AgTr  3 (CTrjTn  + Dj,)”1  c'(Q)AeTr3 

where  m is  the  standard  deviation  of  the  prediction  and  C(Ag)  is  the  variance 
(or  mean  square  value)  of  the  anomaly  block  being  estimated.  Since  we  are 
dealing  with  mean  anomalies,  all  covariances  dealing  with  Ag  values  have 
been  computed  by  numerical  integration  of  the  point  covariance  functions  over 
the  mean  anomaly  block. 

For  this  initial  test, predictions  were  made  for  a 5°  equal  area  block, 
and  a 2°x  2°  block.  A grid  of  8 assumed  data  points  surrounding  the  blocks 
to  be  predicted  was  established  as  is  shown  in  Figure  3.1.  The  predic- 
tions were  carried  out  for  two  assumed  observation  heights  (250  km  and 
850  km)  and  for  several  accuracy  specifications  on  the  ST/3r  values. 

The  predicted  accuracies,  with  respect  to  an  ellipsoidal  reference 
model,  is  given  in  Table  3.1  for  the  5°  block  and  Table  3.2  for  the 
2°  x 2°  block. 


Table  3. 1 

Accuracy  of  the  Predicted  5°  Anomaly  Block 
Using  Known  Data  as  Shown  in  Figure  3. 1 




Data  Acc 

— 1 

Observation  Height 
250  km  850  km 

0.0  mgals 
0.1 
1.0 
2.0 

3.6  mgals 
3.6 
4.0 
4.8 

3.1  mgals 
5.3 

8.1 
10.2 

1 

Table  3. 2 

Accuracy  of  the  Predicted  2°x  2°  Anomaly  Block 
Using  Known  Data  as  Shown  in  Figure  3. 1 


Observation  Height 

Data  Acc 

250  km 

850  km 

0.0  mgals 

6„  3 mgals 

8.9  mgals  j 

0.1 

6.4 

14.6 

1.0 

9.6 

17.8 

to 
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Figure  3. 1:  Location  of  5°  and  2°  Anomaly  Blocks  To  Be 

BT 

Predicted  From  Given  - — Values  at  Indicated  Points  (•) 
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In  judging  the  accuracy  of  figures  given  in  Tables  3. 1 and  3.  2,  it  is 
helpful  to  note  that  the  root  mean  5°  anomaly  is  about  18.9  mgals  while 
the  corresponding  value  for  the  2° x 2C  anomaly  is  25.2  mgals. 

The  values  given  in  these  tables  clearly  show  the  improvement  in 
the  anomaly  determination  as  the  data  accuracy  is  improved  and/or  the 
observation  height  is  reduced. 

The  promising  results  shown  in  these  tables  indicate  that  a more 
complete  simulation  study  should  be  performed  trying  to  formulate  a 
procedure  that  closely  relates  to  a real  world  situation.  Such  a sim- 
ulation is  described  in  the  following  three  sections. 


4.  Simulation  of  Data 

We  consider  the  sum  range  rate,  R8  (Martin,  1972),  from  a 
ground  tracking  station  to  the  geosynchronous  relay  satellite  ATS-6 
tracking  the  close  satellite  GEOS-3  . The  R,  observa- 

tions were  simulated,  considering  that  the  earth's  gravity  field  is 
fully  described  by  the  potential  coefficients  of  the  Goddard  Earth 
Model  8 (GEM  8)  (Wagner  et  al,  1975),  complete  to  degree  and  order 
25  along  with  some  resonance  coefficients  to  (30,  28).  This  will  be 
called  the  high  degree  or  the  (30,  28)  P.C.  field,  and  the  R8 
observations  as  R„(30,  28).  The  7508.0  version  of  the  GEODYN 
program  (Martin  et  al,  1975)  was  used  for  generating  these 
observations. 

The  reference  gravity  field  for  describing  the  earth's  normal 
potential  was  taken  as  the  potential  coefficients,  complete  to  degree 
and  order  12,  from  the  GEM  8 field.  This  will  be  called  the  low 
degree  or  the  (12,  12)  P.C.  field,  and  the  Rs  values  generated  in  this 
low  degree  field  will  be  termed  as  R,(12,  12).  The  residual  range 
; rate  Rrc  between  the  relay  and  close  satellites  was  obtained  as: 

(4.1)  Rre  = R,(30,  28)  - R,(12,  12) 

considering  that  the  relay  satellite  is  not  perturbed  over  short 
periods  by  the  potential  coefficients  higher  than  the  low  degree  field. 

The  residual  acceleration  Rrs  between  the  relay  and  close 
satellites  was  obtained  by  the  numerical  differentiation  of  Rr!.  values. 
This  was  done  by  fitting  a natural  cubic  interpolatory  spline  to  Rrc 
values.  The  details  of  the  procedure  will  be  described  later  in 
Section  4. 3. 
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In  judging  the  accuracy  of  figures  given  in  Tables  3. 1 and  3. 2,  it  is 
helpful  to  note  that  the  root  mean  5°  anomaly  is  about  18.9  mgals  while 
the  corresponding  value  for  the  2°x  2°  anomaly  is  25.2  mgals. 

The  values  given  in  these  tables  clearly  show  the  improvement  in 
the  anomaly  determination  as  the  data  accuracy  is  improved  and/or  the 
observation  height  is  reduced. 

The  promising  results  shown  in  these  tables  indicate  that  a more 
complete  simulation  study  should  be  performed  trying  to  formulate  a 
procedure  that  closely  relates  to  a real  world  situation.  Such  a sim- 
ulation is  described  in  the  following  three  sections. 


4.  Simulation  of  Data 

We  consider  the  sum  range  rate,  Re  (Martin,  1972),  from  a 
ground  tracking  station  to  the  geosynchronous  relay  satellite  ATS-6 
tracking  the  close  satellite  GEOS-3  . The  Rs  observa- 

tions were  simulated,  considering  that  the  earth’s  gravity  field  is 
fully  described  by  the  potential  coefficients  of  the  Goddard  Earth 
Model  8 (GEM  8)  (Wagner  et  al,  1975),  complete  to  degree  and  order 
25  along  with  some  resonance  coefficients  to  (30,  28).  This  will  be 
called  the  high  degree  or  the  (30,  28)  P.  C.  field,  and  the  Ra 
observations  as  Rs(30,  28).  The  7508.0  version  of  the  GEODYN 
program  (Martin  et  al,  1975)  was  used  for  generating  these 
observations. 

The  reference  gravity  field  for  describing  the  earth's  normal 
potential  was  taken  as  the  potential  coefficients,  complete  to  degree 
and  order  12,  from  the  GEM  8 field.  This  will  be  called  the  low 
degree  or  the  (12,  12)  P.C.  field,  and  the  Re  values  generated  in  this 
low  degree  field  will  be  termed  as  R„(12,  12).  The  residual  range 
rate  Rrc  between  the  relay  and  close  satellites  was  obtained  as: 

(4.1)  Rro  = R,(30,  28)  - R,(12,  12) 

considering  that  the  relay  satellite  is  not  perturbed  over  short 
periods  by  the  potential  coefficients  higher  than  the  low  degree  field. 

The  residual  acceleration  Rro  between  the  relay  and  close 
satellites  was  obtained  by  the  numerical  differentiation  of  Rrc  values. 
This  was  done  by  fitting  a natural  cubic  interpolatory  spline  to  Rre 
values.  The  details  of  the  procedure  will  be  described  later  in 
Section  4.  3. 


The  inertial  position  and  velocity  vectors,  and  the  latitude,  longitude 
and  height  above  the  ellipsoid,  of  the  relay  and  close  satellites  were  also 
generated  in  the  low  degree  field  using  the  Geodyn  program.  If  at  time  t, 
the  inertial  position  coordinates  of  the  relay  and  close  satellites  are  given 
r » Yr , zr,  Xe,  Yc,  Za  respectively,  then  the  position  vectors 
£r  * £o  the  relay  and  close  satellites,  and  the  respective  unit  position 
vectors  er  and  ea , as  well  as  the  position  vector  rra  from  the  relay  to 
close  satellite,  may  be  computed  from; 


(4.2)  er  = £/  |rr  I = (Xr  i_  + Yr  j_  + Zr  k)  / (X2  + Yr2  + Zr2)Va 

£c  = rc/  |ra  | = (Xa  i_  + Ya  j_  + Zak)  / (X?  + Y?  + Z2)7a 
era=  rrc/  | rr0 1 = ^(Xa  - Xr)  i_  + (Y0  - Yr ) j_  + (Zo-Zr)  / ^(Xo  -X,f+(%- Yf+(Ze- Zrf 


with  usual  notations  for  the  magnitude  of  a vector,  and  the  unit  coordinate 
vectors. 


The  angle  /3  between  the  radial  direction  at  the  close  satellite,  and 
the  direction  from  the  close  satellite  to  the  relay  satellite  may  then  be 
computed  from: 


(4.3)  cos0  = (-  ea)  • (-erc)  = e 0 • ero 

= (xa(Xe-  X,)  +%(Ya-Yr)  +Za(Za-Zr))/  I5  | k.1 


The  anomalous  potential  T with  respect  to  the  low  degree  field  is 
then, using  approximation  (I), related  to  Rra  values  through  the  relation: 

(4*4)  (dT/3r)a  = RPC  / cos  0 

The  residual  gravity  anomalies  Ag'  referred  to  the  low  degree 
field  may  then  be  predicted  as  discussed  in  Section  3.  The  data  required 
for  these  computations  are  the  latitude  <pa , longitude  \o , height  ha  of 
the  close  satellite  along  with  (3  T /S  r)8  at  the  close  satellite  at 
various  times. 


Further  details  about  the  generation  of  this  data  are  described  in 
this  section. 
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4. 1 Residual  10°  and  5°  Equal  Area  Mean  Anomalies 


The  potential  coefficients  in  the  GEM  8 field  are  referred  to  the 
following  values  of  the  gravitational  constant  times  the  mass  of  the  earth  kM, 
the  semi-major  axis  a,  of  the  ellipsoid,  and  the  fully  normalized  potential 
coefficient  of  degree  2 order  0,  C20  as: 

(4.5)  kM  = 3- 986,008  x 1014  m3/sec2 

a,  = 6, 378, 145  m 
Cs,  o = -484-  164,  57  x 10"6 

The  flattening  f of  the  ellipsoid  may  be  computed  (Hajela,  1974  b,  page  7) 
by  iteration  from: 


(4.6) 


f 


= f3_ 
L 2 


J.  + 


co 


2„3 


il± 


2k  M 


( 9 

l't 


\ 


li 

49 


r + 


/ V'-i 


where  to  is  the  rotational  velocity  of  the  earth,  and  the  dynamical  form 
factor  JP  is  given  by: 


Jg  - - C3  , c 

resulting  in  the  value  of  flattening  as: 
(4.7)  f = 1/298.  257 


The  mean  gravity  anomaly  AgPC  of  a block,  bounded  by  geocentric 
latitudes  <&',  , Os'  and  longitudes  XE  , , as  implied  by  the  potential 

coefficients  up  to  degree  N„  4 „ is  given,  after  neglecting  the  zero  degree 
value  of  the  gravity  anomaly,  by  (Desrochers,  1971,  page  13): 


(4.8)  AgPC 


A X 

k M / r2  V 
AA(sino/-  sino/) 


,„-i,  (*-) 


AX  C*. 


1 Os 


P., 


N 

(sino')  coso'  do'  + j — < CN  ,m  (sinmXE  - sinmAw)  - 
SN  fM  (cosmXE  - cosmAw)  J • j PN  (sino')  coso'  do  j , 

Os 


where  AX  = A*  - Art;  and  r is  the  geocentric  radius  to  the  center  of  anomaly 
block.  CN  fM  and  Sm  >m  are  the  fully  normalized  potential  coefficients  of 
degree  n and  order  m,  the  Ctf.o  being  reduced  by  Cs,0and  C4(6;  and Pn,m (sin<p') 
are  the  fully  normalized  associated  Legendre  functions  of  degree  n and 
order  m. 

The  residual  mean  gravity  anomaly  Ag7  of  a block  referred  to  the 
low  degree  field  was  computed  as: 

(4.9)  Ag'  = Ag(30,  28)  - Ag(12,  12) 

These  gave  the  expected  values  of  the  residual  mean  gravity  anomalies, 
against  which  the  predicted  values  may  be  compared.  10°  equal  area 
and  5C  equal  area  residual  mean  anomalies  were  considered  according 
to  the  subdivision  scheme  for  equal  area  blocks  described  in  (Hajela, 

1974  a,  page  28).  These  will  be  referred  to  as  simply  10°  and  5° 
anomalies,  and  it  will  be  understood  that  we  are  discussing  the  re- 
sidual anomalies  referred  to  the  (12,  12)  P.C.  field. 

Figures  4. 1 and  4. 2 show  the  anomaly  blocks  which  will  be 
considered  in  this  paper.  The  anomaly  blocks  have  been  numbered  for 
easy  reference,  separately  for  10°  and  5°  blocks.  The  particulars  of 
these  anomalies  are  given  in  Tables  4.1  and  4.2. 


Table  4. 1 

10°  Equal  Area  Mean  Residual  Anomalies 


1 

No. 

AS  i A? 

Ag(30, 28) 
mgals 

4?(12,12) 

mgals 

Ag' 

mgals 

! i 

50 

40 

263 

277 

- 6.82 

- 2.89 

-3.93 

2 

40 

30 

252 

264 

3.24 

5.33 

-2.09 

3 

40 

30 

264 

276 

- 5.39 

0.04 

-5.42 

4 

40 

30 

276 

288 

- 7.24 

-14.60 

7.36 

! 5 

30 

20 

251 

262 

6.93 

2.13 

4.80  j 

1 6 

30 

20 

262 

273 

- 2.09 

5.44 

-7.53 

7 

30 

20 

273 

284 

- 2.67 

- 8.83 

6.15 

8 

20 

10 

257 

267 

14.22 

11.01 

3.20  I 

9 

1 

20 

L_^_J 

267 

278 

18.64 

15.92 

2.72 

Table  4.2 

5°  Equal  Area  Mean  Residual  Anomalies 


<P° 

\° 

A.W 

X? 

Ag(30, 28) 
mgals 

Ag(12, 12) 
mgals 

Ag' 

mgals 

1 

35 

258 

264 

- 0.13 

7.11 

- 7.24 

2 

35 

264 

270 

- 7.03 

3.98 

-11.01 

3 

35 

252 

258 

5.86 

2.82 

3.04 

4 

35 

258 

264 

- 1.25 

5.47 

- 6.72 

5 

35 

264 

-11.25 

2.82 

-14.07 

6 

35 

270 

276 

1.24 

- 4.23 

5.47 

7 

■ : 

251 

257 

3.74 

4.22 

8 

25 

257 

262 

5.79 

1.72 

9 

25 

262 

268 

- 7.93 

5.19 

-13.13 

fl 

25 

268 

273 

- 6.24 

- 6.58 

11 

25 

■ . 

257 

262 

15.79 

6.56 

9.23 

12 

25 

262 

268 

4.97 

11.04 

- 6.06 

4.2  Residual  Range  Rate  Between  Relay  and  Close  Satellites 

The  area  of  investigation  for  the  prediction  of  10°  and  5°  anomalies 
was  centered  on  30°  N and  270°  E and  was  a portion  of  the  area  discussed 
in  (Hajela,  1974  a).  The  choice  of  the  initial  epoch  for  defining  the  inertial 
coordinate  system,  and  the  computation  of  the  inertial  position  and  velocity 
vectors  for  the  ATS-6  and  GEOS-3  satellites  at  the  starting  time  for  each 
arc  was  described  in  that  reference,  and  is  therefore  not  being  repeated. 
The  location  of  the  arcs  with  reference  to  the  10°  and  5°  anomaly  blocks 
have  been  shown  in  Figures  4.1  and  4.2  respectively.  The  direction  of 
movement,  ascending  or  descending,  of  the  subsatellite  point  of  the  close 
satellite  has  also  been  indicated.  The  longitudinal  spacing  was  the  same 
for  the  ascending  arcs  and  the  descending  arcs.  In  each  case,  it  was 
about  6° , which  was  roughly  one-half  of  the  size  of  the  10°  equal  area  block. 

The  particulars  of  the  arcs  are  given  in  Table  4. 3. 
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Table  4.3 

Satellite  Arcs  Used  For  Prediction  of  10°  & 5C  Residual  Anomalies 


Arc 

No. 

Ascending/ 

Descending 

Starting  Time* 

Duration  of 
Arc 

First  & Last  Subsatellite  Pts. 
For  Close  Satellite 

Arc 

Day 

Hour 

Min. 

Sec. 

Min. 

Sec. 

<0c 

X? 

(Pc 

x? 

1 

0 

09 

04 

00 

22 

20 

57 

306 

- 9 

250 

2 

# 

0 

10 

44 

30 

19 

20 

61 

292 

7 

233 

3 

\ 

0 

18 

44 

30 

22 

20 

- 8 

299 

59 

240 

4 

0 

20 

25 

00 

19 

30 

-14 

277 

47 

236 

5 

* 

4 

08 

44 

00 

20 

10 

50 

305 

-12 

261 

6 

4 

10 

23 

30 

21 

30 

57 

293 

- 6 

239 

7 

4 

18 

26 

00 

18 

50 

7 

304 

61 

247 

8 

> 

4 

20 

03 

30 

20 

30 

- 9 

287 

54 

239 

9 

* 

/ 

2 

08 

53 

50 

21 

20 

54 

306 

-10 

256 

10 

* 

2 

10 

33 

00 

19 

20 

61 

299 

7 

240 

11 

2 

18 

34 

40 

21 

20 

- 3 

303 

60 

242 

12 

V 

2 

20 

14 

30 

19 

40 

-11 

282 

50 

238 

13 

6 

10 

10 

50 

23 

30 

60 

306 

- 9 

244 

14 

"v 

6 

19 

51 

00 

22 

20 



-13 

295 

56 

241 

* in  elapsed  time  from  the  initial  epoch  defining  the  inertial  coordinate  system. 


The  pattern  of  arcs  shown  in  Figure  4.1  would  be  available  in  a period 
of  6 days,  if  2 ascending  and  2 descending  arcs  are  observed  each  day.  The 
sum  range  rate  (Rs)  observations  were  simulated  for  about  20  minutes  for 
each  arc  till  the  subsatellite  point  for  the  close  satellite  traversed  the  area 
roughly  bounded  by  latitudes  10°  S.  to  60c  N.  and  longitudes  240°  E.  to  300°  E. 
The  Rs  observations  were  simulated  on  a tape  separately  in  the  high  degree 
(30,  28)  P.  C.  field  and  the  low  degree  (12,  12)  P.C.  reference  field.  This 
was  done  by  using  the  Geodyn  program  in  the  data  reduction  mode,  with 
input  as  zero  Rs  observations  every  10  seconds  in  the  Geodyn  binary 
format  (Martin  et  al,  1975,  page  C 26.16),  for  the  duration  of  arc  as  shown 
in  Table  4.3.  The  residuals,  with  their  signs  reversed,  then  simulated 
Rs  (30,  28)  and  Rs  (12,  12)  observations  at  time  interval  of  10  seconds, 
depending  upon  whether  the  high  degree  or  low  degree  field  was  used  in  the 
Geodyn  program.  The  residual  range  rate  Rr0  between  the  relay  and  close 
satellites  was  obtained  as  in  equation  (4.1): 

Rr a = Rs  (30,  28)  - Rs  (12,  12) 
and  written  out  on  another  tape  in  separate  files  for  each  arc. 
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The  geodetic  coordinates  of  the  ground  tracking  station  for  simulating 
Rs  observations  for  all  arcs  were  assumed  as: 


<p  = 35°  11 ' 56*"096 
X = 277°  07  ' 27*"899 
h = 840*  52  meters 


giving  a hypothetical  tracking  station  in  the  vicinity  of  Rosman,  North  Carolina. 


4.3  Residual  Acceleration  Between  Relay  and  Close  Satellites 


We  first  determine  a natural  cubic  interpolatory  spline  fitting  a finite 
number  of  Rr0  values,  and  then  compute  the  first  derivative  of  the  spline 
at  any  of  the  above  points,  or  at  an  interpolated  point  between  them.  The 
fitted  spline  would  differ  depending  on  the  number  of  Rre  point  values  used, 
eg.  10,  15,  or  20  points,  and  whether  we  use  these  points  at  every  10,  20, 
or  30  seconds  apart.  The  value  of  R^  would  then  also  differ  depending  on 
the  number  of  Rre  values  used,  and  their  time  interval. 


If  we  consider  n values  of  Rr( 


R1,Rs,...,Rn  at  time  1 1 , t a , ...,  tn. 


then  we  define  the  spline  S(t),  continuous  along  with  its  first^(S  (t)  ) and 
second  (S##(t)  ) derivatives  on  the  closed  interval ) Rx,  Rr.  j » and  fulfilling 
the  conditions: 


(4.11) 


(a)  S(tj)  - R i ; i = 1,  2,  . .. , n 

(b)  S"(tx)  = S"  (tn)  =0 

(c)  S(t)  is  a cubic  polynomial  on  each  interval 


l”t  * , 1 1 j ; i — If  2,  • • • f n-1 


The  equation  of  the  spline  under  the  above  conditions  is  given  by: 
(For  further  details,  see  (Shampine  and  Allen,  1973,  page  54)  . The 
algorithm  described  below  is  based  on  that  reference.) 


(4.12) 


RTO  ■ s-(‘)  = 557  ^ (t't,) 


S'/+1  At,  N« 


) (t  - t ,) 


(JBj.  - 

\ At , 


s !'nt , 


) (tui  - t ) ; 


i lj  2f  • • • f n-1, 


i 


i 

I 

: 


! 


1 

I 

: i 

;-i 


where: 

(4.13)  tj  st^tlfl  and  Att=tul-t! 


We  are  interested  in  the  first  derivative  of  the  spline,  which  is  obtained 
by  differentiating  equation  (1.12),  and  is  given  by: 


(4.14) 


R(t)  = S/(t)  = - (t1+1-t)2  + 


2 At , 


(t-tt)£ 


+ 


Pt»i~  Ri 
At , 


i 1,  2,  ...,  n-1 


The  second  derivatives  (S")  of  the  spline  needed  in  equation  (4.14) 
are  computed  by  the  recursive  relation: 


(4.15)  S,"  = PjS/'  + t,  ; i=2,  ....  n 

starting  from  Sn7/  = 0 ; after  we  have  computed  the  coefficients  pl  and  t4  by 
the  following  recursive  relations  (Shampine  and  Allen,  1973,  page  57); 
starting  with  pB  = ts  = 0 in  view  of  equations  (4. 11  b)  and  (4. 15) : 

(4.16)  P,.,  = -1/(^4-  p,  4 2(l  + ) 


(4.17) 


■ t ■)/(f^p’+2(i+£f0); 

i = 2,  . . . , n ; and  where 


(4.18) 


d: 


G 

At 


(3  ~ Ri  _ R 1 - R1-1  ^ 

1 V Att  At , / 


The  variation  of  it(t)  in  equation  (4.14)  was  tested  using  R,.c  values 
in  2 arcs  (Arcs  1 and  8 in  Figure  4.1),  when  the  spline  was  determined 
usings,  10,  15,  20,  25  It,.  c values,  which  maybe  10,  20,  30,  40,  50,  60 
seconds  apart.  The  computation  of  It  (t) , i.e.  ftp  c , was  in  each  case  done 
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at  points  10  seconds  apart,  and  the  n points  (n  = 5,  10,  ....  25)  were  equally 
distributed  on  each  side  of  the  computation  point.  The  last  condition  of  the 
distribution  of  points  for  the  spline  was  of  course  limited  by  the  first  and 
last  Rre  value  in  the  arc,  when  the  computation  points  were  in  the  beginning 
or  end  of  the  arc. 

It  was  found  that  if  the  spline  was  fitted  on  Rrs  values  at  20,  30,  40 
or  60  seconds  apart  instead  of  Rrc  values  at  every  10  seconds  apart,  there 
was  a noticeable  change  in  the  computed  R,.0  values.  This  may  be  seen  in 
Table  4.4  where  the  computation  points  are  51  and  55,  i.e.  8 min.  20  sec. 
and  9 min.  respectively  from  the  beginning  of  the  arc;  number  of  points 
(Rr0 ) used  to  define  the  spline  are  20  and  15,  and  the  time  interval  of 
these  R,.s  points  are  10,  20,  30,  40  and  60  seconds. 

Table  4. 4 

Variation  Of  Residual  Acceleration  Due  To  Time  Interval  Of 
Points  For  Fitting  The  Cubic  Spline 


# Pts.  In 
Spline 

Time 
Interval 
(Sec. ) 

Rr  s W 

mgals 

Compn.  Pt.  In  Arc  1 

Compn.  Pt.  In  Arc  8 

51 

55 

51  55 

20 

10 

0.54559 

0.19493 

-0.43189  -1.09578 

20 

20 

0. 54624 

0.19324 

-0.43234  -1.09459 

20 

30 

0.54614 

0.19352 

-0.43314  -1.09415 

20 

40 

0.54640 

0.19355 

-0.43229  J -1.09507 

20 

60 

0.54311 

0. 19255 

-0.43696  -1.09215 

15 

10 

0.54560 

0.19490 

-0.43184  -1.09581 

15 

20 

0.54607 

0.19449 

-0.43214  -1.09544 

15 

30 

0.54601 

0. 19380 

-0.43126  -1.09535  | 

15 

40 

0.54314 

0.19416 

-0.43115  -1.09518 

15 

L 

60 

0.53866 

0.19988 



-0.43733  , -1.08928 

• 

It  is  therefore  advisable  to  use  the  Rr.  values  at  10  seconds  apart 
to  determine  the  spline.  The  change  in  Rre  values  with  this  time  interval, 
but  with  different  values  of  n = 5,  10,  15,  18,  20,  25  is  shown  in  Table  4.  5. 
We  note  that  there  is  little  change  in  Rr0  after  n exceeds  15. 
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Table  4. 5 

Variation  of  Residual  Acceleration  Due  to  No.  of  Points 
Used  For  Fitting  The  Cubic  Spline 


1 

# Pts.  In 
Spline 

Time 
Interval 
(Sec. ) 

RPc  in 

mgals 

1 

Compn.  Pt. 

In  Arc  1 

Compn.  Pt 

. In  Arc  8 

51 

55 

51 

55 

i 

25 

10 

0.54559 

0.19494 

-0.43189 

-1.09578 

20 

10 

0.54559 

0.19493 

-0.43189 

-1.09578 

18 

10 

0. 54559 

0.19495 

-0.43188 

-1.09578 

15 

10 

0.54560 

0.19490 

-0.43184 

-1.09581 

10 

10 

0.54573 

0.19721 

-0.43100 

-1.09324 

5 

10 

0.52741 

0.22564 

-0.48189 

-1.05332 

The  computation  of  Rr0  for  all  arcs  was  finally  done  after  fitting  the 
cubic  spline  to  18  Rro  values  at  10  seconds  time  interval,  i. e.  for  a total 
time  span  of  3 minutes.  With  this  fitted  spline,  Rr,  was  computed  at  10 
seconds  time  interval  for  6 points  in  the  central  1 minute.  The  Rr=  values 
were  thus  computed  for  one  minute  time  span  on  each  occasion  using  Rr0 
values  for  3 minutes,  centrally  straddling  the  one  minute  period  over  which 
Rr0  was  being  computed.  Only  in  the  beginning  first  minute  and  the  last 
one  minute,  or  the  balance  portion  of  the  minute,  the  3 minute  Rrc  values 
did  not  encompass  the  computation  points  for  Rrc  centrally.  It  was 
however  found  that  due  to  condition  (4.11b),  the  Rrc  values  in  the  beginning 
and  the  end  of  the  arc  are  much  less  sensitive  to  the  number  of  points 
used  for  fitting  the  spline. 

The  Rre  values  were  written  out  on  tape,  one  file  for  each  arc,  at 
10  seconds  time  interval  throughout  the  duration  of  arc  as  given  in  Table  4.  3. 

4.4  Radial  Derivative  of  the  Disturbing  Potential 

The  ephemeris  for  the  relay  and  close  satellites  were  written  out  on 
tape  in  the  low  degree  (12,  12)  P.C.  field  by  the  Geodyn  program  as 
described  in  the  beginning  of  Section  4.  The  ephemeris  provided  the 
inertial  position  and  velocity  vectors,  and  the  geodetic  latitude  , longitude 
and  height  of  both  satellites  at  10  seconds  time  interval.  The  information 
for  all  arcs  being  processed  in  a run  was  available  on  the  same  file  on 
the  tape,  along  with  headings  and  arc  summaries  suitable  for  printer 
output.  The  ephemeris  information  was  extracted  out  separately  for  each 
arc,  and  written  on  separate  files  on  another  tape. 
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With  the  inertial  position  coordinates  of  both  relay  and  close  satellites 
available  in  each  arc  at  10  seconds  time  interval,  the  cosine  of  the  angle 
between  the  radial  direction  at  the  close  satellite  and  the  direction  from 
the  close  satellite  to  the  relay  satellite  could  be  computed  at  any  time 
according  to  equation  (4.  3).  Again,  with  the  Rrc  values  available  at  the 
same  times,  as  described  in  the  end  of  Sec.  4.  3,  the  radial  derivative 
of  the  disturbing  potential  (3T/9r)c  at  the  close  satellite  could  be  computed 
using  equation  (4.4). 

The  epliemeris  information  for  each  arc,  and  the  residual  accel- 
eration Rr0  values  for  each  arc,  could  then  be  merged  together  to  provide 
all  the  data  required  for  predicting  the  residual  anomalies.  The  pertinent 
information  consists  of  time;  geodetic  latitude  <p0 , longitude  X5  and 
height  he  of  the  close  satellite;  and  (3T/dr)e  values.  This  was  written 
out  at  10  seconds  time  interval,  one  file  for  each  arc  for  the  duration 
given  in  Table  4.  3,  on  the  final  data  tape  to  be  used  for  the  prediction  of 
residual  anomalies. 

The  algorithms  for  the  prediction  of  residual  anomalies  have  been 
described  in  Sections  2 and  3.  The  actual  prediction  tests  using  the  final 
data  tape  will  be  described  in  Sections  5 and  6. 

5.  Prediction  of  10°  Equal  Area  Residual  Anomalies 

We  will  now  examine  the  prediction  of  10c  residual  anomalies  using 
the  (5T/or)c  values  simulated  in  the  previous  section.  The  expected 
value  of  the  anomalies,  E(Ag'),  has  been  given  in  Table  4.1.  The 
estimated  value,  , of  a particular  anomaly  is  computed  from; 

(5.1)  4?'  = CA6'Tr  (CTrTr  + Df1  T r 

where  Tr  =(3T/3r)0  is  a vector  of  a finite  number  of  (&T/c*r)c  values 
used  for  estimating  the  anomaly,  with  an  assumed  noise  matrix  D and 
the  covariance  matrix  C Tr  Tr  , and  C^g'Tr  represents  the  transpose  of 
the  vector  of  covariances  between  the  anomaly  block  being  estimated 
and  the  Tr  vector. 

The  computation  of  CA*'Tr  and  CTr  Tr  was  discussed  in  Section  2 . 

We  now  assume  a diagonal  noise  matrix  D with  all  elements  the  same, 
and  equal  to  the  square  of  the  assumed  standard  deviation  of  data,  Tr  . 

We  will  later  consider  the  effect  on  the  prediction,  of  various  assumed 
values  of  the  standard  deviation  of  the  data. 
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The  number  of  data  points,  (9T/Sr)c  values  in  the  vector  Tr , 
depends  on  the  spherical  distance.  Ip,  from  the  center  of  the  anomaly 
block  within  which  the  data  is  considered.  Denoting  the  geodetic  latitude 
and  longitude  of  the  center  of  anomaly  block  as  <p'6  and  , and  those  of 
a (ST/9r)0  data  point  as  <pc , Xc  ; the  value  of  ip  is  given  (Heiskanen  and 
Moritz,  equation  2-168): 


(5.2) 


Ip  = cos' 


sin  (fig  sin  <pa 


+ cos  (fig  cos  oB 


The  number  of  data  points  will  also  depend  on  whether  we  consider  each 
data  point  at  10  seconds  time  interval  along  the  arc,  or  use  some  other 
time  interval,  say  30  seconds  or  1 minute.  And,  finally,  the  number  of 
data  points  will  depend  on  the  longitudinal  spacing  of  arcs.  It  may  be 
recalled  from  Figure  4.1  that  the  longitudinal  spacing  of  ascending,  as 
well  as  descending,  arcs  was  about  one-half  of  the  10°  anomaly  block. 


We  will  examine  the  variation  of  the  prediction  depending  on  the 
noise  matrix,  or  the  assumed  standard  deviation  of  the  data;  and  the 
variation  due  to  number  of  data  points  depending  primarily  on  the 
spherical  distance  ip,  and  also  on  the  time  interval,  and  the  longi- 
tudinal spacing  of  arcs. 


The  quality  of  the  collocation  will  be  examined  by  the  standard 
deviation,  b^  , of  the  estimated  anomaly;  as  well  as  the  anomaly 
discrepancy,  e,,-,  and  for  comparison  only,  the  correlation  coefficient, 
p,  of  the  predicted  and  expected  values  of  the  anomalies.  These 
statistics  may  be  computed  from: 

(5.3)  = Co  -C^g'Tr  (CTrTr  +D)  C Tr 

where  C0  is  the  variance  of  the  anomalies  of  that  block  size. 

(5.4)  eAs'  = 4k'  ~ E(Ag') 

(5.5)  p = (z  Ag'  E ( Ag') /n  (( Z A& ' 2/n)V£  (l  (E  ( Ag')  f / n)V®  ) 

where  n is  the  number  of  anomalies  being  predicted. 

A low  root  mean  square  (R.  M.  S. ) value  of  the  anomaly  discrepancy 
as  compared  to  the  mean  standard  deviation  of  the  predicted  anomalies, 
and  a high  value  of  correlation  coefficient  would  prove  the  feasibility  of 
the  prediction  method. 


5. 1 Variation  Due  to  Change  in  Assumed  Standard  Deviation  of  Data 


During  the  simulation  of  data  in  Section  4,  we  did  not  consider  any 
observational  errors  on  the  R,  observations.  Also,  no  uncertainty  was 
ascribed  to  the  potential  coefficients  describing  the  high  degree  or  low 
degree  gravity  field,  or  to  the  epoch  parameters  of  the  arcs.  We  there- 
fore simulated  perfect  data  of  (ST/9r)c  values  at  the  positions 
(*e , Xe  , h- ) of  the  close  satellite  at  various  times,  except  for  any 
inherent  model  inadequacy  and  degeneration  in  the  procedures  of  obtaining 
(3T/Sr)c  values  from  simulated  R,  observations. 

We  may  now  consider  this  by  treating  the  Tj.  data  to  have  an  assumed 
standard  deviation.  6 cases  were  considered  with  the  standard  deviation  as 
0.0,  0.1,  0.5,  1.0,  2.0,  and  5.0  mgals.  We  present  in  Table  5.1  the  pre- 
dicted anomaly,  Ag',  and  the  standard  deviation  of  the  predicted  anomaly, 
on  6 ' , for  two  10°  anomalies  numbered  3 and  6 in  Figure  4. 1.  We  also 
introduce  different  number  of  data  points  by  letting  the  spherical  distance, 
il),  vary  as  5°,  7.5°,  10°,  15°.  The  time  interval  of  data  points  was  kept 
as  1 minute  throughout,  except  for  an  additional  test  with  30  seconds  for 
if,  - 5°.  All  the  14  arcs  in  Figure  4.1  were  considered. 

We  first  note  that  with  assumed  standard  deviation  of  data  as  0.0 
mgals,  though  aA6'  is  obtained  satisfactorily,  the  predicted  anomaly 
Ag'  is  greatly  in  error.  By  comparing  equations  (5. 1)  and  (5.3),  we 
conclude  that  the  simulated  data  of  (5T/9r)c  values  are  not  usable  with 
standard  deviation  as  0.0  mgals.  The  prediction  is  improved  with  standard 
deviation  of  data  as  0.1  mgals,  but  is  still  unstable  for  varying  spherical 
distance  y.  On  the  other  hand,  as  the  standard  deviation  of  the  data  is 
increased  beyond  1 mgal,  the  predicted  anomaly,  Ag , is  damped  off; 
and  this  is  very  clearly  noticable  with  the  standard  deviation  of  data  as 
5 mgals. 


To  examine  this  further,  the  predicted  anomaly  is  plotted  against 
the  assumed  standard  deviation  of  data  in  Figure  5.1,  where  we  present  three 
cases  of  0 — 5°,  with  time  interval  of  data  as  30  seconds  and  1 minute;  and 
ifi=  7.  5°  with  time  interval  of  data  as  1 minute.  The  results  are  plotted  not 
only  for  10°  anomalies  numbered  3 and  G given  in  Table  5. 1 and  whose 
expected  values  were  negative  , but  also  for  anomalies  4 and  7,  whose 
expected  values  were  positive.  The  predicted  anomaly  for  standard 
deviation  of  data  as  0.0  mgals  has  been  omitted. 

We  clearly  notice  in  all  cases  the  damping  of  predicted  anomalies 
as  the  standard  deviation  of  data  is  increased  beyond  1 mgal.  The  some- 
what erratic  behavior  of  the  predicted  anomaly  with  standard  deviation  of 
data  as  0. 1 mgal  is  also  noted  in  some  cases.  The  optimum  standard  de- 
viation of  the  simulated  (3T/dr)a  data  appears  from  Figure  5. 1 to  be 
0.5  mgals.  (This  will  be  found  still  more  noticeable  in  Figure  6.1  for 
for  the  case  of  5°  anomalies.) 
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Table  5. 1 

Prediction  of  10°  Residual  Anomalies  Showing  Variation  Due  To 
Assumed  Standard  Deviation (S.  D.)  of  Data,  Spherical 
Distance  (^°)  From  Center  of  Anomaly  Block, 

And  Time  Interval  of  Data  Points 


Spherical 
Radius  1 1>° 
& Time 
Interval 
of  Data 

Assumed 
S.  D.  of 
Data 
(mgals) 

Ar 

Exp< 

10m.  # 3 
scted  Value 

Anom.  # 6 
Expected  Value 
= -7.5  mgals  1 

# Data 
Pts. 

V**  & 

(mgals) 

# Data 
Pts. 

V & 

(mgals) 

0.0 

3.2 

- 74.2 

3.2 

- 2.6 

Ip 

= 5° 

0.1 

3.6 

- 4.7 

3.6 

-8.3 

0.5 

18 

4.3 

- 7.4 

16 

4.2 

- 8.3 

30  sec. 

1.0 

4.6 

- 7.1 

4.6 

-7.6 

2.0 

5.4 

- 6.0 

5.4 

-6.3 

5.0 

7.0 

- 3.0 

7. 1 

- 2.9 

0.0 

3.6 

- 38.6 

3.5 

- 9.6 

Ip 

= 5° 

0.1 

4.1 

- 10.0 

3.8 

- 8.4 

0.5 

8 

4.4 

- 7.5 

6 

4.4 

- 8.1 

1 

min. 

1.0 

5.0 

- 6.8 

5.0 

- 7.3 

2.0 

6. 1 

- 5.0 

6.2 

- 5.1 

5.0 

7.6 

- 1.7 

7.6 

- 1.6 

0.0 

3.2 

- 12.6 

3.2 

-30.0 

Ip 

11 

•<1 

C7I 

O 

0.1 

3.4 

- 5.9 

3.4 

-13.2 

0.5 

19 

4.2 

- 6.4 

19 

4.1 

-10.0 

1 

min. 

1.0 

4.9 

- 6.2 

4.9 

- 7.3 

2.0 

5.9 

- 5.1 

5.9 

- 4.9 

5.0 

7.4 

- 2.1 

7.4 

- 1.8 

0.0 

3.0 

95.9 

3.1 

-15.7 

Ip 

= 10° 

0.1 

3.3 

- 7.6 

3.3 

-10.6  | 

0.5 

31 

3.9 

- 7.0 

30 

3.9 

- 9.6 

1 

min. 

1.0 

4.6 

6.  5 

4.8 

- 7.5 

2.0 

5.  8 

CM 

LO 

1 

5. 9 

- 5.0 

5.0 

7.4 

- 2.1 

7.4 

- 1.8 

0.0 

2.9 

-134.7 

2.9 

-70.0 

lp 

= 15° 

0.  1 

3.2 

- 6.6 

3.3 

-12.2 

0„  5 

60 

3.8 

- 7.0 

75 

3.8 

- 9.6 

1 

min. 

1.0 

4.4 

- 7.3 

4.5 

-7.5 

2.0 

5.4 

- 6.4 

5.5 

- 5.3 

5‘° 

7., 

- 3.2 

7.1 

-2.2 

$ . = Standard  deviation  (mgals)  of  predicted  1 0"  residual  anomaly 
A^r '=  Predicted  10c  residual  anomaly  (mgals) 


9 
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if/  =5° 

Data  Interval  30  Sec. 


Y Axis  Predicted  ANOM.  ( MGALS) 


X Axis  S.D.  of  Data  (MGALS) 


= 5° 


Data  Interval  1 Min. 


Figure  5. 1:  Variation  of  predicted  10°  residual  anomalies  due  to 
change  in  assumed  standard  deviation  of  data 
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The  prediction  of  other  10°  anomalies  would  thus  be  examined  with 
standard  deviation  of  data  as  0.5  mgals,  even  though  the  agreement  of 
predicted  anomaly  with  the  expected  value  may  appear  to  be  closer  in 
some  cases  with  standard  deviation  of  data  as  1 mgal. 


5. 2 Variation  Due  to  Change  in  Spherical  Distance  From  the  Center  of 


Anomaly  Block 

The  values  for  anc*  Ag*  for  anomalies  3,  4,  6 and  7 are  given  in 
Table  5.2  for  the  spherical  distance  ip  = 5°,  7.5°,  10°  and  15°.  The  time 
interval  of  data  was  1 minute  in  all  cases,  and  the  standard  deviation  of 
the  data  was  taken  as  0.5  mgals. 


Table  5. 2 

Prediction  of  10°  Residual  Anomalies  Showing  Variation  Due  To 
Spherical  Distance  (ip°)  From  Center  of  Anomaly  Block. 
Assumed  Standard  Deviation  of  Data  “0.5  mgals. 

Time  Interval  of  Data  = 1 Minute. 


Anom.  # 3 
E(Ag') 


Anom.  # 4 
E(Ag') 

= 7.4  mgals 


Anom.  # 6 
E(Ag') 

= -7. 5 meals 


Anom.  # 7 
E(Ag') 

= 6. 2 meals 


= Standard  deviation  (mgals)  of  predicted  10°  residual  anomaly 
4?'  = Predicted  10°  residual  anomaly  (mgals) 

We  note  that  the  improvement  in  the  standard  deviation  of  the 
predicted  anomaly  is  very  slight,  from  4.4  mgals  to  3.8  mgals,  as  the 
number  of  data  points  considered  (due  to  increasing  ip)  changes  from  6 to  8 
points  at  lp=  5°  to  63  to  86  points  at  Ip  = 15°. 
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The  estimated  anomaly  also  shows  only  a slight  change  due  to  change 
in  spherical  distance  ip.  For  a standard  deviation  of  about  4 mgals, 
the  range  of  Ag  is  1.1,  0.4,  1.9,  0.7  mgals  for  the  4 anomalies  as  ip 
changes  from  5°  to  15°.  The  small  range  in  the  predicted  values  is  also 
shown  in  Figure  5.2,  where  the  predicted  anomaly  Ag'  is  plotted  against 
spherical  distance  ip. 

It  therefore  appears  that  the  prediction  of  the  anomaly  block  is  mainly 
dependent  on  the  data  points  which  are  close  to  the  center  of  the  block,  i.  e. 
within  Ip  = 5°  for  the  10c  blocks.  Additional  data  points  beyond  ip=  5°  do 
not  contribute  any  significant  information  to  the  predicted  anomaly,  or 
significantly  lower  the  standard  deviation  of  the  predicted  anomaly. 

The  location  of  data  points  with  reference  to  the  anomaly  block  may 
be  seen  in  Figures  5.3  (a)  and  (b)0  These  show  the  data  points  considered 
at  time  interval  of  1 minute  for  the  prediction  of  10°  anomaly  number  6, 
when  Ip  is  taken  as  5°  and  7.5C.  The  number  of  data  points  increase  from 
6 to  19,  but  the  additional  13  points  lie  on  the  periphery  or  in  the  comers 
of  the  block  and  do  not  contribute  any  significant  information  as  seen  from 
Table  5.2. 

Further  tests  for  the  prediction  of  10°  residual  anomalies  need  then 
be  done  only  with  data  points  within  a spherical  distance  of  5°  from  the 
center  of  the  block. 


5.  3 Variation  Due  to  Change  in  Time  Interval  of  Data  Points 

Figures  5.  3 (c)  and  (d)  show  the  location  of  data  points  for  the  pre- 
diction of  10°  anomaly  number  6,  each  for  ip  = 5C,  but  with  time  interval 
of  data  points  as  1 minute  and  30  seconds  respectively.  The  number  of 
data  points  increases  from  6 to  16.  We  present  in  Table  5.3  the  standard 
deviation  of  the  predicted  anomaly  and  the  predicted  anomaly  Ag' 
due  to  change  in  time  interval  of  data  points  for  1 minute  to  30  seconds. 
The  information  pertains  to  10°  anomalies  numbered  3,  4,  6 and  7. 
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PREDICTED  ANOMALY  (MGALS) 


Figure  5.2:  Variation  of  predicted  10°  residual  anomalies  due  to 
change  in  spherical  distance  (t°)  from  center  of  anomaly 
block  up  to  which  data  was  considered.  Data  interval 
was  1 minute.  Standard  deviation  of  data  was 
assumed  as  0.  5 mgals. 


Table  5. 3 

Prediction  of  10°  Residual  Anomalies  Showing  Variation  Due  To 
Time  Interval  of  Data  Points.  Spherical  Distance  From 
Center  of  Anomaly  Block  = 5°.  Assumed  Standard 
Deviation  of  Data  = 0.5  mgals. 


Time 
Interval 
of  Data 
Pts. 
(sec.) 

Anom.  # 3 
E(Ag') 

= -5. 4 mgals 

Anom.  # 4 
E(Ag') 

= 7.4  mgals 

Anom.  # 6 
E(Ag') 

= -7.  5 mgals 

Anom.  # 7 
E (Ag  ) 

= 6.  2 mgals 

# 

Pts. 

A 

V 

Ag 

# 

Pts. 

A 

A*  * 

Ag 

# 

Pts. 

oQ> 

Ag 

# 

Pts. 

A 

4' 

60 

30 

7 

19 

4.4 

4.3 

-7.5 

-7.4 

8 

18 

4.4 

4.2 

8.5 

8.6 

6 

16 

4.4 

4.2 

-8.1 

-8.3 

7 

17 

— 

4.4 

4.3 

6.4 

6.3 

= Standard  deviation  (mgals)  of  predicted  10°  residual  anomaly 
= Predicted  10°  residual  anomaly  (mgals) 


It  is  obvious  that  the  time  interval  of  1 minute  for  the  data  points  is 
adequate  for  the  prediction  of  10°  residual  anomalies. 


5.4  Summary  of  Results 

We  now  present  in  Table  5.4  the  predicted  values,  Ag',  of  9 10° 
residual  anomalies,  and  the  standard  deviation,  , of  the  predicted 
anomalies.  The  anomaly  discrepancy,  eAg',  i.  e.  the  discrepancy 
between  the  predicted  and  the  expected  value  according  to  equation  (5.4)  is 
also  given  in  Table  5.4.  The  expected  value  of  the  anomalies  was  given  in 
Table  4. 1,  and  the  location  of  the  anomalies  in  relation  to  the  satellite 
arcs  was  shown  in  Figure  4. 1.  We  recall  that  the  longitudinal  spacing  of 
ascending,  and  descending  arcs,  was  roughly  one-half  of  the  10°  equal 
area  blocks,  and  that  the  time  interval  of  (dT/dr)c  data  points  was  1 
minute.  We  consider  data  points  within  a spherical  distance  of  5°  from 
the  center  of  anomaly  block,  and  we  assume  the  standard  deviation  of 
data  as  0.  5 mgals. 
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Table  5.  4 

Prediction  of  10°  Residual  Anomalies. 

Spherical  Distance  of  Data  Points  From  Center  of  Anomaly  Block  = 5° 
Time  Interval  of  Data  Points  = 1 Minute.  Assumed  Standard 
Deviation  of  Data  = 0.5  mgals. 


Anom. 

No. 

d? 

<Ps 

SO 

E (Ag') 
(mgals) 

# Data 
Pts. 

A 

V 

(mgals) 

Ag 

(mgals) 

eA*' 

(mgals) 

1 

50 

40 

263 

277 

-3.93 

11 

4.32 

-0.78 

3.15 

2 

40 

30 

252 

264 

-2.09 

8 

4.38 

-0.22 

1.87 

3 

40 

30 

264 

276 

-5.42 

7 

4.41 

-7.48 

-2.06 

4 

40 

30 

276 

288 

7.36 

8 

4.38 

8.46 

1.10 

5 

30 

20 

251 

262 

4.80 

6 

4.34 

0.23 

-4.57 

6 

30 

20 

262 

273 

-7.53 

6 

4.37 

-8.10 

-0.57 

7 

30 

20 

273 

284 

6.15 

7 

4.43 

6.44 

0.29 

8 

20 

10 

257 

267 

3.20 

6 

4.42 

3.41 

0.21 

9 

20 

10 

267 

278 

2.72 

7 



4.26 

4.64 

1.92 

Limits  of  Anomaly  Block  <p°  , North  & South  Latitude  in  Degrees 

X°rt  , X°e  West  & East  Longitude  in  Degrees 
E (Ag" ) Expected  Value  of  10c  Residual  Anomaly 
! u,,'  Standard  Deviation  of  Predicted  10°  Residual  Anomaly 

Ag*  Predicted  10°  Residual  Anomaly 

Anomaly  Discrepancy  = Predicted  Ag^  - Expected  E(Ag‘>) 

The  mean  standard  deviation  of  the  predicted  anomalies  was  4 . 37 
mgals,  while  the  root  mean  square  (R.  M.  S.)  value  of  anomaly  discrepancy 
was  2.21  mgals.  Further,  as  the  R.  M.  S.  value  of  the  expected  anomalies 
was  5.15  mgals,  the  R.  M.  S.  value  of  anomaly  discrepancy  is  acceptable. 
The  correlation  coefficient,  p,  between  the  predicted  and  expected  values 
of  the  9 anomalies  considered,  as  computed  from  equation  (5.  5),  was 
0.915.  We  may  therefore  conclude  that  the  procedures  developed  for  the 
prediction  of  10°  residual  anomalies  are  workable. 

i 

•i  ] 

6,  Prediction  of  5°  Equal  Area  Residual  Anomalies 

We  first  test  out  the  conclusions  arrived  at  in  Sections  5. 1 to  5.  3 
for  10°  residual  anomalies,  regarding  their  applicability  in  the  case  of 
5°  residual  anomalies.  We  then  use  the  considerations,  so  developed, 
regarding  the  standard  deviation  of  the  (cT/Sr)c  data  and  the  number  of 
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data  points  to  be  used,  in  predicting  12  5°  residual  anomalies  shown  in 
Figure  4.2.  A comparison  of  the  predicted  anomalies  against  the  expected 
values  given  in  Table  4.2  would  enable  us  to  examine  the  anomaly  discrep- 
ancies according  to  equation  (5.4),  and  the  correlation  coefficient,  p, 
between  the  predicted  and  expected  values  according  to  equation  (5.  5). 

A low  root  mean  square (R.  M.  S. ) value  of  anomaly  discrepancy  in 
comparison  to  mean  standard  deviation  of  the  predicted  anomalies,  and 
a high  value  of  the  correlation  coefficient,  p,  would  reinforce  our  con- 
clusions of  the  feasibility  of  the  prediction  method  developed  in  this  paper. 


6. 1 Variation  Due  to  Change  in  Assumed  Standard  Deviation  of  Data 

We  again  examine  6 values  of  the  assumed  standard  deviation  of 
(3T/Sr)c  data.  These  were  0.  0,  0.1,  0.  5,  1.0,  2.0  and  5.0  mga^s  as 
in  Section  5. 1.  The  standard  deviation  of  the  predicted  anomaly,  aaS' , 
and  the  predicted  anomaly,  Ag  , are  given  in  Table  6.1  for  two  5° 
residual  anomalies  numbered  5 and  9.  Because  of  the  smaller  area  of 
anomaly  block,  we  consider  the  time  interval  of  data  points  as  30  seconds, 
instead  of  1 minute  in  the  case  of  10°  anomalies.  The  spherical  distance, 

0 , from  the  center  of  anomaly  block  within  which  time  data  points  were 
considered,  was  |/)=2.5°,  3.5°,  5°,  7.5°,  10°. 

We  first  note  from  Table  6. 1 that  for  standard  deviation  of  data  as 
0.0  mgals,  though  is  obtained  satisfactorily,  it  is  net  possible  to  use 
the  data  for  prediction  of  5°  anomalies.  The  somewhat  erratic  behavior 
of  the  predicted  anomaly  with  standard  deviation  of  0. 1 mgals  is  also 
noticeable.  This,  and  the  damping  of  the  predicted  anomaly  for  standard 
deviation  of  data  exceeding  1 mgal,  is  seen  very  clearly  in  Figure  6. 1, 
where  the  predicted  anomaly  has  been  plotted  against  the  assumed 
standard  deviation  of  the  data  for  the  spherical  distance  = 2.5°,  3.5° 
and  5°. 

We  therefore  conclude,  as  in  Section  5.1,  that  for  the  simulated 
(9T/3r  )c  data  in  Section  4,  it  is  appropriate  to  use  the  standard  devi- 
ation of  0.  5 mgals. 


6.2  Variation  Due  to  Change  in  Spherical  Distance  From  the  Center  of 
Anomaly  Block 

The  standard  deviation  of  the  predicted  anomaly,  b-it'  , and  the 
predicted  anomaly,  Ag*,  have  been  shown  in  Table  6.2  for  4 5°  residual 
anomalies  numbered  2,  5,  9 and  12  for  the  spherical  distance  = 2.5°, 
3.5°,  5°,  7.5°,  10°.  The  time  interval  of  data  was  30  seconds  in  all 
cases,  and  the  data  was  assumed  to  have  a standard  deviation  of  0.  5 
mgals. 
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Table  6. 1 

Prediction  of  5°  Residual  Anomalies  Showing  Variation  Due  To 
Assumed  Standard  Deviation  (S.  D. ) of  Data,  and  the  Spherical 
Distance  (0°)  From  Center  of  Anomaly  Block. 


Spherical  Assumed  Anom.  # 5 Anom.  # 9 

Radius  0 S.  D.  of  Expected  Value  = -14. 1 mgals  Expected  Value  = -13. 1 mgals 
Time  Data  

Interval  (mgals)  # Data  A'e  # Data 

of  Data  Pts.  i (mga 


214.0 

- 3.2 

- 13.3 

- 11.6 

- 7.7 

- 2.3 

-242.3 

- 6.0 

- 13.2 

- 13.2 

- 10.4 

- 3.9 


-303.9 

- 9.9 

- 13.1 

- 13.6 

- 11.9 

- 5.9 


0= 

2.5° 

30 

Sec. 

0 = 

3.5° 

30 

Sec. 

0 = 

5° 

30 

Sec. 

0 = 

:7.5° 

30 

Sec. 

0 = 

10° 

30 

Sec. 

Standard  Deviation  (mgals)  of  Predicted  5 Residual  Anomaly 
Predicted  5°  Residual  Anomaly  (mgals) 


Y Axis  Predicted  ANOM.(MGALS)  X Axis  S.D.  of  Data  ( MGALS ) 


Figure  6. 1:  Variation  of  predicted  5°  residual  anomalies  due  to 
change  in  assumed  standard  deviation  of  data 
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Table  6.  2 

Prediction  of  5°  Residual  Anomalies  Showing  Variation  Due  To 
Spherical  Distance  (^°)  From  Center  of  Anomaly  Block. 
Assumed  Standard  Deviation  of  Data  = 0.5  mgals. 

Time  Interval  of  Data  = 30  Seconds. 


— 

Anom.  # 2 
E (Ag') 

= -11.  0 mgals 

Anom.  * 5 
E (Ag') 

= -14. 1 mgals 

Anom.  # 9 
E (Ag') 

= -13. 1 mgals 

Anom.  14  12 
E (Ag') 

= -6. 1 mgals 

# 

Pts. 

V 

aA  ' 

Ag 

# 

Pts. 

V' 

Ag' 

# 

Pts. 

V 

Ag' 

# 

Pts. 

A 

(T  5' 
a* 

— 

Ag' 

2.5° 

4 

11.2 

-10.6 

4 

10.9 

-13.3 

5 

10.  6 

-18.2 

ii.  i 

- 3.6 

3.5° 

9 

11. 1 

- 9.7 

8 

10.9 

-13.2 

6 

10.6 

-18.1 

8 

10.7 

- 5.7 

5.0° 

20 

10.9 

- 8.0 

18 

10.7 

-13.1 

15 

10.3 

-17.8 

17 

10.5 

- 5.9 

7.5° 

44 

9.8 

- 9.8 

40 

9.6 

-12.0 

40 

9.4 

-21.6 

33 

9.6 

- 7.1 

10° 

1 

80 

9.7 

-11.1 

74 

9.  5 

. 

-12.1 

69 



9.3 

-21.4 

68 

9.3 

-10.  1 

= Standard  Deviation  (mgals)  of  Predicted  5°  Residual  Anomaly 


= Predicted  5°  Residual  Anomaly  (mgals) 


There  does  not  appear  to  be  much  difference  in  the  predicted  anomaly, 
or  much  reduction  in  the  standard  deviation  as  the  spherical  distance  is 
increased  from  2.5°  to  10°,  except  for  the  predicted  anomaly  for  b = 2.5° 
for  anomaly  12,  which  may  be  explained  by  only  2 data  points  being  con- 
sidered. Because  of  the  likelihood  of  insufficient  data  points,  we  may  not 
consider  $ = 2.5°,  but  \b  = 3.5°  and  5.0C  appear  to  be  equally  suitable  for 
predicting  the  5°  residual  anomalies.  Further  tests  will  accordingly  be 
reported  with  i/j=  3.5°.  However,  when  all  the  12  5 residual  anomalies 
were  finally  predicted  considering  spherical  distances  $=  3. 5 , 5°,  7.5""; 
the  latter  sets  were  found  to  be  slightly  better.  This  will  be  further 
elaborated  in  Section  6.4. 


6.3  Variation  Due  to  Change  in  Time  Interval  of  Data  Points 

We  present  the  values  of  by  and  Ag'  in  Table  6.3  for  5°  residual 
anomalies  numbered  2,  5,  9 and  12  with  time  interval  of  data  points  as  30 
and  10  seconds.  Data  points  were  considered  within  the  spherical  distance 
of  3.  5°  and  the  standard  deviation  of  data  was  assumed  as  0.  5 mgals. 


-44- 


Table  6.3 

Prediction  of  5°  Residual  Anomalies  Showing  Variation  Due  To 
Time  Interval  of  Data  Points  Spherical  Distance  From 
Center  of  Anomaly  Block  = 3.  5°. 

Assumed  Standard  Deviation 
of  Data  = 0.5  mgals. 


Time 

Interval 

of 

Data 

Pts. 

(Sec.) 

Anom.  # 2 
E (Ag') 

= -11.  0 mgals 

Anom.  # 5 
E (Ag') 

= -14.1  mgals 

Anom.  # 9 
E (Ag') 

= -13. 1 mgals 

Anom.  # 12 
E (Ag') 

= -6. 1 mgals 

# 

Pt& 

V 

Ag' 

# 

Pts. 

V 

A A * 

Ag 

# 

Pts. 

V 

aa  * 

Ag 

# 

Pts. 

6 

.A  ^ 

Ag 

30 

10 

9 

29 

11. 1 
10.9 

- 9.7 
-10.3 

8 

28 

10.9 

10.5 

-13.2 
-11. 1 

6 

22 

10.6 

10.3 

-18. 1 
-16.  8 

8 

28 

. 

10.7 

10.4 

-5.7 

-4.7 

a^'  = Standard  Deviation  (mgals)  of  Predicted  5°  Residual  Anomaly 


Ag'  = Predicted  5°  Residual  Anomaly  (mgals) 

It  is  obvious  that  the  time  interval  of  30  seconds  is  adequate  for  the 
data  points  for  the  prediction  of  5C  residual  anomalies,  and  no  significant 
advantage  is  gained  by  the  increase  in  the  number  of  data  points  due  to  the 
reduction  in  the  time  interval  of  data  points  from  30  seconds  to  10  seconds. 
The  location  of  data  points  has  been  shown  in  Figures  6.2  (a)  and  (b)  for 
residual  anomaly  number  9 with  the  time  interval  of  data  as  30  and  10 
seconds  respectively. 

We  recall  from  Figure  4.2  that  the  longitudinal  spacing  of  arcs 
was  roughly  the  same  as  the  5°  anomaly  block  size,  while  in  the  case  of 
prediction  of  10°  anomalies,  the  longitudinal  spacing  of  arcs  was  one-half 
of  the  size  of  10°  block.  Additional  satellite  arcs  were  then  simulated 
till  the  longitudinal  spacing  of  ascending,  and  descending,  arcs  was  one- 
half  of  the  5°  anomaly  block.  The  location  of  data  points  with  this  increased 
density  of  arcs  has  been  shown  in  Figure  6.2  (c)  for  anomaly  number  9 
with  0=  3.5°  and  time  interval  of  data  points  as  30  seconds.  In  Figure 
6.2  (d)  we  consider  the  same  case  but  with  time  interval  of  data  points 
as  60  seconds. 

The  values  of  a ■ ^ and  Ag'for  anomaly  number  9 have  been  shown 
in  Table  6.4  for  the  2 cases  in  Figures  6.2  (c)  and  (d> , i.e.  for  longi- 
tudinal spacing  of  arcs  one-half  of  the  5°  block  size  with  time  interval  of 
data  as  30  seconds  and  60  seconds.  The  values  for  the  case  in  Figure 


6.2  (a),  i.e.  spacing  of  arcs  same  as  the  5 block  size  with  time  interval 
of  data  as  30  seconds  have  also  been  shown  for  comparison. 


r 


262°  E. 


268°E. 


262°E. 


268°  E 


30°N. 


25°N. 


Arc  Long.  Spacing  = Block  Size 
(a) Data  Interval  30  Seconds 


Arc  Long.  Spacing  = Block  Size 
( b)  Data  Interval  10  Seconds 


» , 

1=  i 


262°E  268°E. 


262°E.  268°E. 


\ 30° N. 


25°N. 


Arc  Long  Spacing  = | Block  Size 
(c)  Data  Interval  30  Seconds 


Arc  Long  Spacing  = ^ Block  Size 
(d)  Data  Interval  60  Seconds 


Figure  6.2;  Location  of  data  points  for  prediction  of  5°  residual 
anomaly  no.  9.  Spherical  distance  from  center  of 
anomaly  block  = 3°.  5 
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'Fable  6.4 

Prediction  of  5°  Residual  Anomaly  No.  9 Showing  Variation  Due  To 
Change  in  Number  of  Data  Points  With  Increased  Density  of  Arcs. 
Spherical  Distance  From  Center  of  Anomaly  Block  = 3.5°. 
Assumed  Standard  Deviation  of  Data  =0.5  mgals 


Anom.  #9  # Data 

E (Ag')  = -13.1  mgals  Pts. 



A 

°A6' 

( mgals) 

Ag' 

( mgals) 

| 

Longitudinal  Spacing  of  Arcs  = 1/2  of  5°  Block  Size 

1 

Time  Interval  of  Data  Points  = 30  Sec.  (Fig.  6.2(c)) 

16 

10.5 

-17. 1 

Time  interval  of  Data  Points  = 60  Sec.  (Fig.  6.2(d)) 
Longitudinal  Spacing  of  Arcs  = Same  as  5C  Block  Size 

6 

10.6 

-16.  8 

Time  Interval  of  Data  Points  = 30  Sec.  (Fig.  6.2(a)) 

1 

6 

10.6 

-18.1 

From  a comparison  of  Figures  6.2  (a)  and  (c),  it  appears  that  a greater 
density  of  observations  within  the  comparatively  small  size  of  a 5C  anomaly 
block  does  not  significantly  alter  the  predicted  values.  Figure  6.2  (d)  appears 
to  be  basically  comparable  to  Figure  6.2  (a)  as  though  the  number  of  arcs 
have  been  increased,  the  time  interval  was  reduced  proportionately. 


6.4  Summary  of  Results 

We  present  in  Table  6.  5 the  standard  deviation  of  the  predicted  anomaly, 
& , the  predicted  anomaly,  Ag',  and  the  anomaly  discrepancy,  , 
according  to  equation  (5.4)  of  all  the  12  5°  residual  anomalies  shown  in 
Figure  4.2.  The  expected  value  of  these  anomalies  was  given  in  Table  4.2. 
Data  points  were  used  up  to  a spherical  distance  of  3.  5°  with  time  interval 
of  data  points  as  30  seconds.  The  standard  deviation  of  the  data  was 
assumed  as  0.5  mgals. 


Table  6.  5 

Prediction  Of  A Local  Set  Of  5°  Residual  Anomalies. 

Spherical  Distance  Of  Data  Points  From  Center  Of  Anomaly  Block  = 3.5° 
Time  Interval  Of  Data  Points  = 30  Sec.  Assumed  Standard  Deviation 

Of  Data  =0.5  mgals. 


Anom. 

No. 

<Ps 

— 

Aw 

A? 

— 

E (Ag') 
(mgals) 

# Data 
Pts. 

A 

<V 

(mgals) 

aA  * 

Ag 

(mgals) 

(mgals) 

1 

40 

35 

258 

264 

- 7.24 

11 

11.1 

- 1.70 

5.  54 

2 

40 

35 

264 

270 

-11.01 

9 

11.1 

- 9.69 

1.32 

3 

35 

30 

252 

258 

3.04 

6 

10.8 

2.59 

-0.45 

4 

35 

30 

258 

264 

- 6.72 

6 

10.9 

-12. 55 

-5.83 

5 

35 

30 

264 

270 

-14.07 

8 

10.9 

-13.22 

0.85 

6 

35 

30 

270 

276 

5.47 

8 

11.0 

1.57 

-3.  90 

7 

30 

25 

251 

257 

4.  22 

7 

10.6 

1.71 

-2.  51 

8 

30 

25 

257 

262 

1.  72 

7 

11.2 

- 6.96 

-8.68 

9 

30 

25 

262 

268 

-13. 13 

6 

10.6 

-18.08 

-4.95 

10 

30 

25 

268 

273 

- 6.58 

6 

11.4 

- 1.15 

5.  43 

11 

25 

20 

257 

262 

9.23 

8 

11.4 

- 0.28 

-9.  51 

12 

25 

20 

262 

268 

- 6.06 

8 

10.  7 

- 5.  66 

0.40 

Limits  of  Anomaly  Block  , <2°  North  & South  Latitude  in  Degrees 

Xw  > Xe  West  &•  East  Longitude  in  Degrees 
E (Ag  ) Expected  Value  of  5°  Residual  Anomaly 

Standard  Deviation  of  Predicted  5C  Residual  Anomaly 

Ag  Predicted  5C  Residual  Anomaly 

Anomaly  Discrepancy  = Predicted  Ag  - Expected  E(Ag  ) 

The  mean  standard  deviation  of  the  predicted  anomalies  was  11.0 
mgals,  while  the  root  mean  square  (RMS)  anomaly  discrepancy  was  5. 1 
mgals.  The  latter  value  was  therefore  acceptable,  considering  that  the 
R.  M.  S.  value  of  expected  anomalies  was  8.2  mgals.  The  correlation 
coefficient  of  the  predicted  and  the  expected  values  of  the  12  anomalies, 
as  computed  from  equation  (5.5),  was  0.81.  The  anomaly  discrepancy 
and  the  correlation  coefficient  indicate  satisfactory  prediction  of  5° 
residual  anomalies.  However,  we  wish  to  examine  if  these  could  be 
improved  further. 

The  12  5°  residual  anomalies  were  predicted  again  using  data  up 
to  spherical  distances  $ = 5°  and  7.5°.  The  time  interval  of  data,  and 
the1  assumed  standard  deviation  of  the  data  was  retained  the  same  as  in 
the  case  of  Table  6.  5 with  b = 3.  5 . The  results  of  these  tests  are  given 
in  Table  G.6,  showing  the  anomaly  discrepancies,  and  the  standard 

deviations,  cl6'  , of  the  12  predicted  anomalies. 
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Table  6.  6 

Prediction  of  a Local  Set  of  5°  Residual  Anomalies  Showing  Variation 
of  Spherical  Distance  (ip)  of  Data  Points  From  Center  of 
Anomaly  Block.  Time  Interval  of  Data  Points  = 30  Sec. 
Assumed  Standard  Deviation  of  Data  =0.5  mgals. 


Anom. 

No. 

>P  = 3.5° 

iP=5° 

lp  = 7.5° 

# Data 
Pt& 

A 

oy 

(mgals) 

( mgals) 

# Data 
Pts. 

(mgals) 

€u«' 

Oil  gals) 

# Data 
Pts. 

V 

(mgals) 

(mgals) 

1 

11 

11.1 

5.  5 

21 

10.8 

4.3 

46 

9.8 

1.2 

2 

9 

11.1 

1.3 

20 

10.9 

3.0 

44 

9.8 

1.2 

3 

6 

10.8 

-0.4 

16 

10.7 

-1.0 

42 

9.5 

-1.0 

4 

6 

10.9 

-5.8 

16 

10.7 

-4.6 

42 

9.6 

-5.3 

5 

8 

10.9 

0.8 

18 

10.7 

1.0 

40 

9.6 

2.1 

6 

8 

11.0 

-3.9 

18 

10.8 

-3.3 

40 

9.  7 

-1.5 

7 

7 

10.6 

-2.5 

17 

10.3 

-2.0 

34 

9.5 

-0.3 

8 

7 

11.2 

-8.7 

16 

10.8 

-5.3 

40 

10.0 

-3.0 

9 

6 

10.6 

-5.0 

15 

10.3 

-4.7 

40 

9.4 

-8.4 

10 

6 

11.4 

5.4 

16 

11.2 

5.2 

38 

10.1 

5.6 

11 

8 

11.4 

-9.5 

18 

11.2 

-8.0 

32 

10.4 

-3.9 

12 

8 

10.7 

0.4 

17 

10.5 

0.2 

33 

9.6 

-1.1 

1 

cx^g'  = Standard  Deviation  of  Predicted  5°  Residual  Anomaly 

= Anomaly  Discrepancy  = Predicted  Anomaly  - Expected  Value  of  Anomaly 


The  improvement  in  standard  deviation  with  increase  in  the  spherical 
distance,  ip  , is  only  slight.  The  anomaly  discrepancy  even  shows  worsening 
in  some  cases,  e.g. , anomalies  numbered  3,  5,  9,  10,  12;  and  only  a 
marginal  improvement  in  anomalies  numbered  2 and  4.  However,  there  is 
an  improvement  when  we  consider  the  prediction  of  the  whole  set  of  12 
anomalies,  as  the  spherical  distance  of  data  points  increases  from  3. 15° 
to  7.5°.  This  is  seen  more  clearly  in  Table  6.7. 
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Table  6.  7 

Prediction  of  a Local  Set  of  12  5°  Residual  Anomalies  (Figure  4.2)  - 
Summary  of  Results.  Time  Interval  of  Data  Points  = 30  Sec. 
Assumed  Standard  Deviation  of  Data  = 0.5  mgals. 


Spherical 

Distance 

Standard  Deviation  of 
Predicted  Anomalies 
in  mgals 

Predicted  - Expected  Value 
of  Anomalies 
in  mgals 

Correlation 

Coefficient 

Min. 

Max. 

Mean 

Min. 

Max. 

R.  M.  S. 

(Eqn.  (5.  5)) 

3.5 

RSI 

11.4 

mm 

-9.5 

5.5 

0.814 

5.0 

m 

11.2 

ES39 

-8.0 

5.2 

4.14 

0.869 

7.5 

9.4 

10.4 

9.75 

-8.4 

5.6 

3.73 

0.913 

Hence,  though  the  prediction  of  individual  5°  residual  anomalies 
does  not  show  significant  improvement  in  using  data  points  beyond  a spherical 
distance  of  3.  5°  from  the  center  of  anomaly  block,  it  would  be  a better 
strategy  while  predicting  a local  set  of  5°  residual  anomalies  to  use  data 
points  up  to  a spherical  distance  of  7. 5°  from  the  center  of  each  anomaly 
block. 


7.  Effect  of  Uncertainty  in  Epoch  Parameters  on  the  Prediction  of  Anomalies 


During  the  simulation  of  sum  range  rate  observations,  Rs , in 
Section  4,  the  inertial  position  and  velocity  vectors  at  the  initial  epoch  of 
each  arc  were  the  same  for  generating  the  Ra  (30,  28)  and  R,  (12,  12) 
observations.  The  residual  range  rate  between  the  relay  and  close 
satellites,  Rr0  , computed  in  equation  (4.1);  and  subsequently  the  residual 
acceleration,  R.0  , and  the  (9T/9  r)  data  points  used  for  the  prediction  of 
10°  and  5°  residual  anomalies  in  Sections  5 and  6,  had  the  implicit 
assumption  that  the  epoch  parameters  for  the  arcs  were  known  without  any 
uncertainty.  This  allowed  the  testing  of  procedures  for  the  prediction  of 
anomalies  without  any  'aliasing'  due  to  uncertainties  in  the  a-priori 
determined  epoch  parameters,  i.e.  the  inertial  position  and  velocity 
vectors  at  the  initial  epoch  of  each  arc.  In  practice,  the  epoch  parameters 
are  determined  from  ground  tracking,  and  by  first  carrying  out  'inner 
iteration'  (Martin  et  al,  1975)  for  each  arc  with  the  available  observations. 
Therefore,  the  epoch  parameters  are  known  only  approximately,  and  we 
would  now  examine  how  this  would  affect  the  prediction  of  anomalies  with 
procedures  tested  in  Sections  5 and  6. 


We  consider  this  by  using  different  epoch  parameters  for  generating 
the  R,  (12,  12)  observations  from  those  used  for  generating  the  R„  (30,  28) 
observations.  We  then  obtain  Rrc  , Rrc  and(dT/dr)c  values  as  outlined 
in  Section  4.  If  we  now  use  these  changed  (BT/c  r)c  data  points  for  the 
prediction  of  10°  and  5°  residual  anomalies,  we  would  notice  the  effect  on 
the  prediction  of  the  uncertainties  in  the  epoch  parameters  by  the  amount 
of  differences  introduced  in  generating  the  R3  (12,  12)  observations.  The 
changes  in  the  epoch  parameters  were  introduced  only  in  one  coordinate 
of  the  inertial  position  and  velocity  vectors  for  each  arc  to  give  one  test. 

Thus  there  were  6 test  cases  in  which  the  effect  on  the  prediction  was 
examined  for  change  in  one  of  the  X0,  Y0,  Z0,  X0,  Y0,  Z0  coordinates 
of  the  close  satellite  at  the  initial  epoch  for  each  arc.  The  epoch  parameters 
of  the  relay  satellite  were  not  changed. 

We  first  present  the  predicted  values  of  10°  and  5°  anomalies  with 
the  changed  (Bt/B  r)0  data,  and  then  examine  what  interpretation  may  be 
given  to  the  predicted  values,  so  obtained. 

7. 1 Effect  on  the  Prediction  of  10°  Residual  Anomalies 


We  investigate  the  effect  on  the  prediction  of  10°  anomaly  number  6. 
The  location  of  the  (^T/3r)c  data  points  used  in  Section  5 for  predicting 
this  anomaly  was  shown  in  Figures  5.  3 (a)  and  (c)  within  the  spherical 
distance  of  5°  from  the  center  of  the  anomaly  block,  and  with  the  time 
interval  of  data  as  1 minute.  This  gave  the  number  of  data  points  as  6, 

2 each  for  Arcs  1,  8 and  13,  as  may  be  seen  from  Figure  4. 1.  We  now 
change  the  position  and  velocity  vector  of  the  close  satellite,  one  coordi- 
nate at  a time,  in  generating  the  Rs  (12,  12)  observations  for  each  of 
Arcs  1,  8 and  13.  After  repeating  the  procedures  in  Section  4 for  these 

3 arcs,  we  will  get  6 changed  data  points  but  at  practically  the  same  loca- 
tions as  in  Figures  5.3  (a)  and  (c),  because  of  only  small  changes  in  the 
epoch  parameters,  as  explained  below. 

We  designate  the  G cases  of  changed  (3T/3r)»  data  by  letters  A to 
F.  The  position  coordinates  were  each  reduced  by  10  meters  in  their 
absolute  value  in  Arcs  1,  8 and  13  giving  cases  A,  B and  C.  Similarly,  the 
velocity  coordinates  were  each  reduced  by  1 cm/sec  in  their  absolute  value 
in  Arcs  1,  8 and  13  giving  cases  D,  E and  F.  The  change  by  10  meters  in 
the  position  coordinates  and  by  1 cm/sec  in  the  velocity  coordinates  was 
considered  to  give  a realistic  representation  of  the  uncertainty  in  the  a- 
prior  knowledge  of  the  epoch  parameters.  Further,  as  the  absolute  values 
of  position  and  velocity  coordinates  were  reduced  but  as  some  of  these 
coordinates  were  positive  and  some  negative,  the  final  effect  was  a random 
increase  or  decrease  of  the  coordinates  in  the  considered  arcs  1,  8 and  13. 
The  increase  is  indicated  by  a + sign,  and  the  decrease  by  a - sign  in 
Table  7. 1 for  each  arc  in  cases  A to  F.  The  (?T/Br)e  values  for  the  6 data 
points  in  Figures  5,  3(a)  and  (c)  are  shown  in  Table  7.2  for  cases  A to  F.  The 
values  used  in  the  original  case  in  Section  5,  when  the  epoch  parameters  for 
the  arcs  1,  8 and  13  were  not  changed,  are  also  shown  in  Table  7.2. 
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Table  7. 1 

Changes  in  Epoch  Parameters  for  Arcs  1,  8 & 13. 


Case 

Designation 

Close 
Satell  ite 
Coordinate 

ARC  1 

ARC  8 

ARC  13 

A 

X0 

+ 10  m 

+ 10  m 

+ 10  m 

B 

Y, 

- 10  m 

+ 10  m 

- 10  m 

C 

?o 

- 10  m 

+ 10  m 

- 10  m 

D 

X0 

- 1 cm/sec 

+ 1 cm/sec 

- 1 cm/sec 

E 

Yo 

- 1 cm/sec 

- 1 cm/sec 

- 1 cm/sec 

F 

z3 

+ 1 cm/sec 

- 1 cm/sec 

+ 1 cm/sec 

Table  7.2 

Changed  Data  Points  Due  to  Change  in  Epoch  Parameters 
of  Close  Satellite  For  Prediction  of  10°  Residual 
Anomaly  No.  6.  Spherical  Distance  of  Data 
Points  From  Center  of  Anomaly  Block = 5°. 

Time  Interval  of  Data  Points  = 1 Minute. 


The  standard  deviation  of  the  predicted  anomaly,  , and  the 
predicted  value,  ^g  , of  10°  residual  anomaly  number  9 are  shown  in 
Table  7.  3 using  the  6 data  points  in  Table  7.3  for  each  of  the  6 cases 
A to  F.  The  values  have  been  tabulated,  assuming  that  the  standard 
deviation  of  the  data  points  may  be  0.0,  0. 1,  0.5,  1.0,  2.0  and  5.0 
mgals.  The  corresponding  values  for  the  original  case  in  Section  5, 
when  the  epoch  parameters  were  not  changed,  were  given  in  Table  5. 1 
for  il>  = 5°  and  time  interval  of  data  points  as  1 minute. 


Table  7.  3 

Prediction  of  10u  Residual  Anomaly  No.  6 Showing  Variation  Due  to 
Change  in  Epoch  Parameters  of  Close  Satellite.  Spherical 
Distance  of  Data  Points  From  Center  of  Anomaly 
Block  = 5°.  Time  Interval  of  Data  Points  = 1 Minute. 

No.  of  Data  Points  = 6. 


All  Values  Are  in  Mgals 

= Standard  Deviation  of  Predicted  10s  Residual  Anomaly  No.  6 


&5  - Predicted  10°  Residual  Anomaly  No.  G.  E (Ag  ) =-—7.5  Mgals 


We  notice  that  the  standard  deviation  of  the  predicted  anomaly  is  the 
same  for  all  cases  A to  F,  and  the  original  case  in  Table  5. 1,  for  any 
particular  assumed  standard  deviation  of  the  data  points.  This  is  so  as 
the  location  of  the  data  points  remains  practically  the  same  in  all  cases 
because  of  only  small  changes  in  the  epoch  parameters,  and  there  is  thus 
no  particular  change  in  any  of  the  quantities  in  equation  (5.3).  We  also 
notice  as  in  Section  5. 1 that  the  value  of  the  predicted  anomaly  is  not  usable 
for  standard  deviation  of  data  as  0.0  mgals.  We  have  to  now  interpret 
what  standard  deviation  may  be  assumed  for  the  data  to  predict  the  10c 
residual  anomaly.  This  will  be  done  in  Section  7.3.  We  may,  however, 
note  that  the  standard  deviation  of  2 mgals  for  the  data  gives,  in  general, 
a close  agreement  of  the  predicted  value  in  cases  A to  F with  the  ex- 
pected value  of  -7.5  mgals. 
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7.2  Effect  on  the  Prediction  of  5°  Residual  Anomalies 


We  investigate  the  effect  of  uncertainty  in  epoch  parameters  of  the 
arcs  in  the  prediction  of  5°  residual  anomaly  number  9.  Figure  6.  2 (a) 
shows  the  location  of  6 (9  T/9  r)c  data  points  used  for  predicting  this 
anomaly  in  Section  6.  There  were  3 data  points  each  from  arcs  8 and 
13,  as  may  be  seen  from  Figure  4. 2,  and  these  were  at  time  interval  of 
30  seconds,  and  within  a spherical  distance  of  3.  5°  from  the  center  of 
anomaly  block. 

We  use  the  same  configuration  of  data  points  but  with  the  epoch 
parameters  of  the  close  satellite  changed  during  generation  of  R,  (12,  12) 
observations,  as  shown  in  Table  7.1  for  arcs  8 and  13  under  cases  A 
to  F.  The  changed  value  of  the  (3T/9  r)0  data  points  is  shown  in  Table  7.4 
for  cases  A to  F along  with  the  original  case  in  Section  6,  when  the  epoch 
parameters  of  arcs  8 and  13  were  not  changed. 


Table  7.4 

Changed  Data  Points  Due  to  Change  in  Epoch  Parameters  of 
Close  Satellite  For  Prediction  of  5°  Residual  Anomaly  No.  9. 
Spherical  Distance  of  Data  Points  From  Center 
of  Anomaly  Block  = 3.5°.  Time  Interval 
of  Data  Points  = 30  Seconds. 


Arc 

No. 

Time  in 
Arc 

(5  T/9  r)e  in  mgals 

Min. 

Sec. 

Original 
Case 
(Sec.  6) 

Case 

A 

Case 

B 

Case 

C 

Case 

D 

Case 

E 

Case 

F 

8 

11 

Ka 

3.06 

3.59 

0.17 

2.05 

0.86 

2.53 

8 

11 

Kg 

3.35 

3.78 

0.36 

2.34 

1. 12 

2.89 

8 

12 

m 

3.42 

3.  76 

0.34 

2.41 

1.  18 

7'  ■ 

13 

11 

S3 

2.75 

4.90 

5.58 

3.80 

13 

11 

rj 

2.58 

4.  86 

5.41 

3.77 

13 

12 

20 

1.85 

2.17 

4.58 

5.01 

I 


The  5°  residual  anomaly  number  9 was  predicted  using  the  changed 
data  points  shown  in  Table  7.4  under  cases  A to  F,  and  the  values  of  <v 
and  are  given  in  Table  7.5.  The  standard  deviation  of  data  was 
assumed  to  be  0.0,  0.1,  0.5,  1.0,  2.0  and  5.0  mgals.  The  values  of 
a * and  Ag*  for  the  original  case  in  Section  6,  when  the  epoch  param- 
eters of  arcs  8 and  13  were  not  changed,  were  given  in  Table  6. 1 for 
ip  = 3.5°  and  time  interval  of  data  points  as  30  seconds. 


Table  7.5 

Prediction  of  5°  Residual  Anomaly  No.  9 Showing  Variation  Due  To 
Change  in  Epoch  Parameters  of  Close  Satellite.  Spherical 
Distance  of  Data  Points  From  Center  of  Anomaly 
Block  = 3. 5°.  Time  Interval  of  Data 
Points  = 30  Sec.  No.  of 
Data  Points  = 6. 


I 


Assumed 

Case  A 

Case  G 

Case  C 

Case  D 

Case  E 

Case  F 

S.  D.  of 

L 

Data 

V 

V 

*4' 

V 

V 

A*  * 

El 

& 

0.0 

7.7 

-45.5 

7.7 

197.3 

7.7 

677.6 

7.7 

-38.9 

7.7 

361.6 

7.7 

0.1 

9.8 

-26.7 

9.8 

- 8.6 

9.8 

55.8 

9.8 

-19.9 

9.8 

22.8 

9.8 

0.5 

10.6 

-25.7 

10.6 

-36.3 

10.6 

-20.0 

10.6 

10.6 

-18.4 

10.6 

- 15.1 

1.0 

11.0 

-23.$ 

11.0 

-33.9 

11.0 

-20.6 

11.0 

11.0 

-18.1 

11.0 

- 13.0 

2.0 

12.1 

MfJlM 

12.1 

-25.0 

12.1 

-15.4 

12.1 

12.1 

-13.5 

12.1 

- 9.6 

5.0 

13.8 

Ba 

13.8 

- 8.8 

13.8 

- 5.4 

13.8 

Bufl 

13.8 

- 4.7 

13.8 

- 3.4 

All  Values  Are  in  Mgals 

V * Standard  Deviation  of  Predicted  5°  Residual  Anomaly  No.  9 
- Predicted  5°  Residual  Anomaly  No.  9.  E(Ag')  - -13.1  Mgala 


We  notice  as  in  Section  7. 1 that  the  standard  deviation  of  the  predicted 
anomaly  is  the  same  in  all  cases  A to  F,  and  the  original  case  in  Table  6.1, 
for  a particular  assumed  standard  deviation  of  the  data,  as  the  location  of 
the  data  points  remains  practically  the  same  due  to  only  small  changes  in 
the  epoch  parameters  of  the  arcs.  Further,  the  predicted  anomalies  are  not 
usable  for  standard  deviation  of  data  as  0.0  mgals.  We  also  note,  as  in 
Section  7. 1,  that  there  is  generally  a close  agreement  of  the  predicted 
anomaly  with  the  expected  value  of  -13. 1 mgals,  when  the  standard  deviation 
of  the  data  is  assumed  as  2 mgals. 

7.3  Interpretation  of  Results 

The  predicted  values  of  10°  residual  anomaly  number  6 and  5°  residual 
anomaly  number  9 are  shown  in  Table  7.6,  after  extracting  these  values  for 
the  standard  deviation  of  data  as  2 mgals  from  Tables  7. 3 and  7.  5. 


Table  7.6 

Predicted  Residual  10°  Anomaly  No.  6 & 5°  Anomaly  No.  9 
Showing  Variation  Due  to  Change  in  Epoch  Parameters 
of  Close  Satellite.  Assumed  Standard  Deviation 
of  Data  = 2 Mgals.  Predicted  Anomaly  Values 
Extracted  From  Tables  7.2  & 7.5 


Residual 
Anomaly  No. 

Predicted  Anomaly  ir  Mgals 

S.  D.  of 

Predicted  Anom. 
(Mgals) 

Case  A 

Case  B 

Case  C 

Case  D 

Case  E 

Case  F 

10°  Anom.  # 6 
E (Ag')  = -7.5  Mgals 

m 

-14.4 

- 8.4 

11 

m 

m 

For  Cases  A to  F 
6.2 

5°  Anom.  # 9 
E (Ag')  = -13.1  Mgals 

-17.4 

-25.0 

-15.4  1 

-11.9 

-13.5 

mm 

For  Cases  A to  F 
12.1 

We  find  from  Table  7.  6 that  both  for  the  10  ' anomaly  number  6 and  the 
5°  anomaly  number  9,  the  predicted  anomalies  in  cases  A to  F lie  almost 
within  one  standard  deviation  (last  column  of  Table  7.6)  of  the  expected  value 
(first  column  of  Table  7.  6).  It  therefore  appears  that  even  when  there  is 
uncertainty  up  to  10  meters  and  1 cm/sec.  respectively  in  the  inertial 
position  and  velocity  coordinates  at  the  initial  epoch  for  the  arcs,  the  pre- 
dicted anomaly  would  be  within  one  standard  deviation  of  the  expected  value. 
The  uncertainty  in  the  epoch  parameters, _ i.e.  the  inertial  position  and 
velocity  coordinates  X0  , Y„,  Z0,  X0,  Y,  , Z0  at  the  initial  epoch,  for 
the  arcs  up  to  10  meters  and  1 cm/sec  respectively,  shows  up  as  spurious 
variations  in  the  (3  T/3  r)0  data  with  a standard  deviation  of  about  2 mgals. 

We  recall  that  in  Sections  5. 1 and  6. 1,  we  found  that  the  simulated  (3T/dr)c 
data  appeared  to  be  impressed  with  a standard  deviation  of  about  0.  5 mgals, 
when  the  epoch  parameters  for  the  arcs  were  implicitly  considered  to  be 
known  without  any  uncertainty.  Now  when  the  epoch  parameters  have 
uncertainties,  but  are  considered  to  be  known  perfectly,  the  uncertainties 
masquerade  as  a larger  standard  deviation  of  the  (3  T/3  r)~  data. 

This  is  further  borne  out  by  examining  the  data  used  by  us  in  Tables 
7.  2 and  7.  4.  The  variation  of  the  data  at  each  of  the  6 data  points,  due  to  the 
uncertainties  in  the  epoch  parameters  of  the  arcs  1,  8 and  13,  has  been 
extracted  from  Table  7.2  and  shown  in  Table  7.7.  The  sample  standard 
deviation  of  the  data  at  each  of  the  6 points  as  computed  from  the  different 
values  in  cases  A to  F has  also  been  shown  in  Table  7.  7. 


Table  7.  7 

Variation  in  (3  T/3  r)c  Data  Due  to  Change  in  Epoch  Parameters. 
Values  Extracted  From  Table  7.2 


The  mean  of  sample  standard  deviations  of  data  computed  from  arcs 
8 and  13  in  Table  7.  1 gave  a similar  value  of  1. 6.  This  value  was  also 
computed  separatoh  for  cases  A to  C,  and  cases  I)  to  F,  in  Table  7.2 
resulting  in  1.7  and  1.4.  It  therefore  appears  reasonable  to  assume  that 
the  uncertainties  of  about  10  meters  and  1 cm/sec  respectively  in  the 
inertial  position  and  velocity  coordinates  of  arcs  at  the  initial  epoch  may 
be  accounted  for  by  assuming  the  standard  deviation  of  (ST/3  r)0  data  as 
2 mgals. 

8.  Conclusions 

The  estimation  of  mean  gravity  anomalies  at  the  surface  of  the  earth 
from  satellite  to  satellite  tracking  sum  range  rate  observations  is  in 
tendency  unstable.  It  involves  (1)  a numerical  differentiation  which  yields 
accelerations  from  the  range  rate  data  and  (2)  a continuation  of  the  radial 
derivative,  T-  , of  the  anomalous  potential,  derived  from  the  accelera- 
tions, downward  to  the  surface  of  the  earth.  Both  procedures  tend  to 
amplify  errors  in  the  observations  and  in  the  model. 

We  propose  a method  that  consists  of  a differentiation  of  an  inter- 
polatary  cubic  spline  derived  from  the  range  rate  data  to  determine 
accelerations  (alternative  methods  are  also  discussed)  and  the  least  squares 
collocation  method  for  the  estimation  of  the  mean  gravity  anomalies  at 
the  surface  of  the  earth  from  the  Tr  - data.  For  a theoretical  unique 
solution  the  gradient,  V_T, , of  the  anomalous  potential  would  have  to  be 
computed  from  the  relative  accelerations,  IVc  , between  the  close  and 
the  relay  satellite.  But  the  determination  of  V Te  from  the  R.c  - values 
would  make  necessary  a simultaneous  tracking  of  the  close  satellite  by 
at  least  three  relay  satellites.  For  the  treatment  of  the  actual  situation 
of  having  only  one  relay  satellite  we  propose  two  approximate  solutions. 

The  approximation  used  in  our  simulation  study  assumed  V Tg  of 
having  only  the  radial  component,  Tr  . The  experiment  chosen  was 
satellite  to  satellite  tracking  between  the  GEOS-3  satellite  (h  ’ 850  km) 
and  the  ATS-G  (h  = 3f>000  km). 

The  simulation  brought  the  following  results: 

— The  proposed  method  (differentiation  of  a cubic  interpolatorv 
spline  followed  by  least  squares  collocation)  is  economical 
and  produced  meaningful  results. 
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— A 10°  residual  anomaly  (with  respect  to  a (12,  12)  reference 

field)  can  be  estimated  witha  standard  deviation  of  about  4.  4 mgals.  The 
input  residual  anomalies  could  be  well  reproduced  within  this 
standard  deviation.  A 5°  residual  anomaly  can  be  estimated 
with  a standard  deviation  of  about  11.0  mgals.  Also  here  the  input 
residual  anomalies  could  be  well  reproduced  within  this 
standard  deviation.  We  are  well  aware,  however,  of  the  fact 
that  the  variance  analysis  and  the  simulation  study  may  not 
give  conclusive  indications  for  the  real  data  processing. 

— The  preliminary  least  squares  collocation  variance  analysis 
shows  that  an  essential  improvement  of  the  above  results  can 
be  achieved  either  by  a dens ificat ion  of  the  data  tracks  or  even 
more  by  using  a close  satellite  with  lower  altitude  (eg.  h = 250  km). 

A reduction  of  the  random  observation  noise  level  does  not  very 
much  contribute  to  an  improvement  of  the  estimation. 

— A very  encouraging  fact  is  that  only  satellite  data  in  a very 
limited  area  are  needed  (inside  a cap  with  a radius  of  about 
7.  5°  spherical  distance  around  the  center  of  the  estimated 
anomaly  block).  This  allows  a very  economical  data  pro- 
cessing. In  addition,  in  contrast  to  least  squares  adjustment, 
least  squares  collocation  allows  the  estimation  of  each  mean 
gravity  anomaly  block  independent  from  one  another,  which 
allows  the  recovery’  of  gravity  anomalies  in  restricted  areas. 

— The  model  error  mainly  caused  by  the  mentioned  approxi- 
mation inherent  in  the  method  could  be  reasonably  controlled 
by  assigning  a standard  deviation  of  ±0.  5 mgal  to  the  noise- 
free  Tr  values. 

— Because  of  their  systematic  influence  orbital  errors  may  cause 
a severe  problem  for  all  satellite  to  satellite  tracking  experi- 
ments. In  our  simulation  errors  of  10  m in  the  position  or 
1 cm/sec  in  the  velocity  were  assigned  to  the  starting 
elements  of  the  close  satellite.  These  errors  might  be 
considered  as  an  upper  limit  of  errors  occuring  in  a real 
experiment  with  dense  ground  tracking  coverage.  The 
estimated  anomalies  were  disturbed  but  their  deviation  from 
the  input  values  is  most  times  still  inside  the  already  shown 
standard  deviation. 

Considering  these  results,  an  application  of  the  presented  method  to 
real  satellite  to  satellite  tracking  observations  is  the  next  logical  step. 
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